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ABSTRACT 


A modified Galerkin's method is developed to analyze the dynamic 
problems of multiple-load-path bearingless rotor blades. The develop- 
ment and selection of functions are quite parallel to CAMRAD procedures. 
This greatly facilitates the implementation of the method into the CAMRAD 
program. A software is developed implementing the modified Galerkin's 
method to determine free vibration characteristics of multiple- load- 
path rotor blades undergoing coupled flapwise bending, chordwise 
bending, twisting, and extensional motions. Results are in the process 
of being obtained by debugging the software. 
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NOMENCLATURE 

Area of cross-section 

Mass centroid offset from elastic axis 

Young 1 s modulus 

Shear modulus 

Reference moment of inertia 

Bending moments of inertia about major and minor 
neutral axes 

Mass polar radius of gyration of blade cross section 
about elastic axis 

Mass radii of gyration about major neutral axis and 
about an axis perpendicular to chord 

Length of the load paths 

Mass per unit length 

Reference mass 

Twisting moment about x-axis 

Bending moments about y and z axes, respectively 

Number of load path 

Number of selected functions 

Radius of the rotor 

Tension 

Elastic displacements in the x,y,z directions, respectively 
Axial force 

Shear forces along y and z directions, respectively 

Mutually perpendicular axis system with 
x - along the blade elastic axis 
y - towards the leading edge 


x,y,z 



VI 


3 Pretwist angle 

p Elastic twist about the elastic axis 

u) Frequencies of vibration 


n 


Blade rotational speed 



1. INTRODUCTION 


The bearingless rotorcraft offers reduced weight, less complexity, 
and superior flying qualities. Bearingless rotor technology was 
successfully applied by Sikorsky in the tail rotors of the Blackhawk 
and S-76 helicopters (Ref. 1). Boeing Vertol built the first successful 
bearingless main rotor (BMR) and flew it in 1978 (Refs. 2 and 3). During 
a recent study for concept definition of an Integrated Technology 
Rotor/Flight Research Rotor (ITR/FRR), thirty-three hub concepts were 
proposed and twenty-one of these were bearingless designs (Ref. 4). 

This suggests that the next generation of rotorcraft will most likely 
be equipped with bearingless rotors. Boeing Vertol is currently design- 
ing and advanced bearingless rotor hub for MBB-105 blade (ABRS Project). 

All practical designs of bearingless rotors include multiple load 
paths and one that was flight tested by Boeing Vertol has three load 
paths. The existence of multiple load paths change the dynamic behavior 
of the helicopter significantly, and also provides a variety of options 
to the designer. Therefore, it is important to develop a capability 
to analyze the helicopters with multiple-load-path bearingless hubs. 
CAMRAD (Ref. 5) is a comprehensive rotorcraft analysis program developed 
at NASA-Ames. This program utilizes recently developed technology to 
analyze the following problems: 

1. Trim solution 

2. Performance, loads and noise 

3. Stability derivatives and handling qualities 
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4. Aeroelastic stability 

5. Vibration and gust response 

CAMRAD does not model multiple load paths and many users need this capa- 
ability in CAMRAD. 

The modification of CAMRAD program to include the multiple load 
path is undertaken under the present program. This modification is 
a major effort requiring several years of concentrated study. The modi- 
fication is undertaken basically in two major phases: 

Phase 1: Extension and validation of CAMRAD procedures 

to include multiple load paths. 

Phase 2: Incorporation of procedures developed under 

phase 1 into the CAMRAD code under consultation 
with Dr. Wayne Johnson. 

The first two years will be spent for Phase 1 studies. As a part 
of this phase, a systematic study is undertaken to extend the CAMRAD 
technologies. CAMRAD uses Galerkin's method (Ref. 5) following 
reference 6. The limitations of this method as exists in CAMRAD are: 

1. Single-load-path modelling can only be done 

2. Axial degree- of- freedom is not included 

in the analysis, and this degree-of -freedom is 
important for multiple- load-path modelling 

3. There is no coupling between bending and 
torsion in the blade modes 

Therefore, a modified Galerkin's method as employed in CAMRAD following 
reference 6 is extended for bearingless multiple- load-path blade. A 
computer program is developed to validate the formulation. The 
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details of the formulation and the listing of the computer program are 
included in this status report. The results are in the process of being 
obtained. 
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2. BASIC EQUATIONS 

The linear, homogeneous, undamped equations of motion for simple 
harmonic free vibrations with frequency to can be written as (Ref. 7) 

-(EAu 1 )' - J2 2 mu - w 2 mu =0 (1) 

-(Tv')' + {(EljSin 2 ^ + EI 2 cos 2 p)v" 

+ (EI 2 -EI 2 )cosgsingw"}" - fi 2 mv - u) 2 mv 
2 2 2 

+ ft mesinfty + oj mesin(5$ + (Q mexsing#)' * 0 (2) 


-(Tw')' + {(EI 1 cos 2 g + EI 2 sin 2 3)w" 

2 2 

+ (EI 2 -EI 2 )cosgsingv"}" - u> mw - to mecosfty 
2 

- (Q mexcosfS#)' = 0 

- (GJd) ' ) ' - oj 2 mk 2 + Q 2 m(k 2 - k 2 )cos_f3<J> 
m 1 

2 

- Q mex( -sinftv' + cosgw' ) 

2 2 2 
+ mesinjiv + u> mesinpv - uj mecosfiw = 0 


rR 


T = 


2 

Q mxdx 


(4) 

(5) 


x 


5 


2 

For algebraic simplicity, the coefficients B^, B^, e A> are assumed 
to be zero in Eqs. (2) to (4). Also, Eq. (1) governing the stretching 
motion is added to introduce the axial coupling associated with the 
multiple-load-path blades. 
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3. BEARINGLESS ROTOR BLADES 


The bearingless main rotor that was flight tested by Boeing Vertol 
has three load paths, viz., two fiberglass flexbeams and a filament 
wound torque tube (Ref. 8). The components of this bearingless main 
rotor are shown in Fig. 1. The outboard-ends of the flexbeams and torque 
tube are connected to the blade through a rigid clevis. The inboard- 
ends of the flexbeams are rigidly connected to the hub. The inboard- 
end of the torque tube is restrained in torsion by the control system 
stiffness and is pinned to the hub by a rod-end bearing (translational 
deflections and bending moments are equal to zero). The idealization 
of such a blade is shown in Fig. 2. Locations 1, 2, 3 in this figure 
correspond to the root of the load paths, clevis and tip of the blade 
respectively. The bearingless main rotor that is currently being 
designed by Boeing-Vertol has two load paths similar to the one shown 
in Fig. 3. The formulation presented here is quite general and the 
computer program is developed for two load paths and can be modified 
easily for increased number of load paths. For the analysis of multiple- 
load-path blades, it is important to establish the equilibrium and the 
compatibility relations across the clevis. 

3.1 Equilibrium Across the Clevis 

Consider the free-body diagram for the clevis as shown in Fig. 4 

and let h and h be the y and z coordinates of the ith load path with 
y . z . 

i i 

reference to a coordinate system located at the blade (Point '0'). 
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DESCRIPTION 


1 

2 

3 

4 

5 
$ 
7 

a 

9 

10 

11 

12 

13 

14 


FLEXURE, INBOARD ATTACHMENT 
FLEXURE 

FLEXURE, OUTBOARD ATTACHMENT 
STEEL HUB 

FLEXURE/HUB ATTACHMENT BOLTS 

TORQUE TUBE, OUTBOARD FITTING 

TORQUE TUBE, OUTBOARD ATTACHMENT 

TORQUE TUBE 

PITCH ARM ASSEMBLY 

PITCH LINK ASSEMBLY 

CLEVIS 

TORQUE TUBE/CLEVIS ATTACHMENT BOLTS 
BLADE/CLEVIS ATTACHMENT BOLTS 
ROTOR BLADE ASSEMBLY 


Fig. 1. Components of a Bearingless Main Rotor 





c UPPER HU3PLA IE r’lT CM ARM n TORQUE TUBE 
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Fig. 3. Bearingless Rotor Concept with Torque Tube Enclosing the Flexibeam 
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Force and moment equilibrium require the following relations to be 
satisfied 


n 



i=l 



M 

x 


l 



h V 
y • z. 

' i i 


h V ) 
z. y . 
i J i 


M 

y 


l 



+ h V 
. x 
i 


M 

z 



i=l 


h V 


x . 

l 


X . 
1 


( 6 ) 


(7) 


( 8 ) 


(9) 


( 10 ) 


(ID 


where n = number of load paths. 


3.2 Compatibility Across the Clevis 

Consider the plane of the clevis as shown in Fig. 5. The 
axial displacement of the ith load path can be written as: 





12 


yields : 


( 12 ) 


consistent with the rigid clevis 


( 13 ) 

( 14 ) 

( 15 ) 

( 16 ) 
( 17 ) 
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4 . MODIFIED GALERKIN'S METHOD 


The basic equations of motion given by Eqs. (1) to (4) can be 
written as: 


du 

dx 

V 

X 

EA 




(18) 

d<fr _ 
dx 

M 

X 

GJ 




(19) 

d 2 v 

- C 12 

M 

y 

i Cu 

M 

z 

(20) 

dx 2 

" EI 1 EI 2 

EI 1 EI 2 

d 2 w 

C 22 

M 

y 

C 12 

M 

z 

(21) 

dx 2 

ei 1 ei 2 

EI 1 EI 2 


- —— n 2 mu = u> 2 mu (22) 


2 2 2 

-r — + 0 m(k -k )cos20<j> 
dx 

0 2 mex( -singv' + cosgw') + ^ 2 mesingv 


2 2 2 2 
u mk $ - a) mesin$v + oj mecosgw 
m 


(23) 


d 2 M 


dx 


Hr - (T 4 ^)' - (ft 2 mcxcosf30) ' 
z ax 


2 2 

= oj mw + ft mecos3<t 


( 24 ) 
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d**M 

- (T 7 ^)' + (ft 2 mexsing<|>) ' 

dx 2 dx 

2 2 2 2 
+ Q mesing<|> - Q. mv = oj mv - oj mesing$ 


(25) 


where 

2 2 

= EI^cos g + EI ? sin g 

C ^2 = (EI ? -EI^)cosgsing 
2 2 

C ^2 = Eltsin g + E^cos g 


A convenient set of functions for Galerkin's method are selected 
with the following properties 

1. Satisfy the root boundary conditions of the load paths 

2. Satisfy the clevis compatibility conditions 

3. Satisfy the blade tip conditions 

The deflections and forces are expanded as finite series in terms of 
Che selected functions as shown below. The choice of functions and 
their modifications are discussed in Chapter 5. 


N 

u = ) a. u 

L ln n 

n=l 


u . 

l 


N 

V 

L 

n=l 


a u 
In n . 
i 


v 


I 


a v 
an n 




n=l 
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N 


v . 

1 


■ l 


a_ v 
2n n. 
1 


n=l 

N 


w 


= l a 3n w n 


n=l 

N 


w . 

l 


a_ w 
3n n . 
i 


n=l 

N 


- I 


a, $ 
4n n 


n=l 


N 


■I 


4n y n . 
i 


n=l 


(26) 


N 


- I 


b. V 
In xn 


n=l 

N 


■ l 


v 

In xn . 
i 


n=l 

N 


■ l 


b 0 M 
2n zn 


n=l 
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z . 

1 


N 

l 

n=l 

N 


b- M 
2n zn . 
1 


■ l 


b, M 
3n yn 


n=l 

N 


-I 


M 

3n yn, 


n*l 

N 


= y b M 
L 4n xn i 


n=l 


N 

y b, M 
/ 4n xn. 
^ i 

n=l 


(27) 


The satisfaction of the following equations is the modified 
Galerkin's method for multiple-load-path blades which is equivalent 
to the quotient method described in Ref. 6. (N = Number of functions 
selected; n = number of load paths) 


I I I J 

j-1 l i=l 


rl 

b,.M . 
3j yj . 
o J 


a_ . (T.w! )' 

3j i J t 
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N 

l 

j-1 


[b-.M" - a_.(Tw.) 

i 3 j 3j 3 


+ a. , (fi"mexcosf5<t . ) ' ] w, dx 
J k 


y (-b_.M . + a T.w! + a . e .xcosB . 4> . ) 

L 3j yj. 3j l j. 4j ii i Y j i 


L- j-1 1 i-1 


- ( -b_ .M + a_ .Tw . + a. .Q mexcos^ . ) w. 

3j y. 3j j 4j k 


J x=l 


’ N 

f n 1 


l 

) (b- .M . - h b. .V . ) - b-.M , 

L 3j yj. z L lj xj. 3j yj 


.1=1 

i-1 



N 

r 

> 

n 

r 

1 

■ l 

- 

l J 

I J 

m i (a 3j w j + a e.cosg.^ )w dx 
o JJ i J J ii 

* 

L j=l 

' i-1 

J 



+ 


2 

uj 



+ a, . ecos34> . )w. dx 
4j J k 


( 28 ) 


where 


k - 1, 2, 3, ... N 
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N n 

l * l 

3-1 i-l 


r^ 


Cb 2j M zj. ‘ a 2j (T i V j. ) ' + a 4j (n Vi xsin ¥j i ) ' 


2 2 
+ a, .S2m.e.sin(3.$. - a_ .£2 m.v. ]v, dx 

4j ii l j . 2 1 l j . k . 

j j j J ii 


f > 


N 

♦ l 

3-1 


rR 


[b,.M . 
2j zj 


+ a, . £2 mesin3$ . 
4 J J 


a_.(Tv! )’ + a. . ( £2^mexsinf3$ . ) ' 
2 j J 4 j J 


a.,!3 mv.] v, dx l 

2j J k ' 


l 

L 3-1 


^ (-b 0 ,M^. + a.,.T.v' - a, .a^m.e.xsing,. ) 


n ,n , ■ , i , v . o . ,«* in , w , r 

2j zj, 2j j j. 4j ii i 3 ± 


v i-l 


(-b-.M 1 . + a_.Tv' - a .£2 mexsing<J).) 
2j zj 2j j 4j j 


J x=i 


N f n 

y J y (b 0 .M . - h b. .V . ) _ b .M . 

L I L 2 j zj. x i lj xj i 2 J zj 


3-1 


i-l 


> v, 


x=i 


= 0 ) 


3-1 


ii 

IJ 


rl 

m.(a_.v. - a . .e . sinfi . 0 . )v, dx 

i 2j j . 4j i i j . k . 
o JJ i J J 11 


i-l 


+ ui 


[J 

L J = 1 


r R 


m(a».v. - a. ,esing<Ji . )v. dx 
2j j 4j Hy j k 


( 29 ) 


where k - 1 , 2, 3, ... N 
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f rl [ , 2 

) ( -b. .M . + £2m.e.x(-a«, . cosg.v' + a_ .sinfJ.w! ) 

L J 0 4 J xj. 11 2j H i j. 3j p i j. 


j-1 L i=l 


2 2 2 2 
+ a_.£2 m.e.sinfj.v. + a, .£2 ra.(k -k )cos2S.A. 
2j ii i j. 4j i m 2 H i*j. 

i i 


r 4 > k dx 
i 


I 

j=l 


rR 


■b, .m'.+ £2 mex(-a- .cos(3v'. + a- .sinpw! ) 
■+J xj 2j J 3j M j 


2 2 2 2 
+ a- .£2 mesinfSv . + a. .£2 m(k -k )cos28<6.)'d>, dx 
2j j 4j m 2 m 1 py j | T k 


N f n 

+< l l [b..M. +h (-b-.M 1 . + a,.T.w! 

L 4 J xj. y. 3j yj 3j i j 

j-1 l i=l 11 1 1 


2 i i 

+ a. .£2 m.e.xcosg.fj). ) - h (-b-.M 1 , + a_.T.v', 

4j ii l J . z . 2j zj . 2j i j . 

l i JJ i JJ i 


- a,. £2 m. e . xsinft . $ . )] - b, .M . 
4J i i P i^J i 4j xj 


' *k > 


x=l 


= 0 ) 


l l L m i (a 4J k i + a 3 j e i C0S P i w j 


Lj-i 


i=l 


i i 


- a_ .e.sing.v. ) <f. dx 
2j i i J . k . 

j i 


+ co 


r f R 2 

L J ," (a 4j k mb + a 3j ecose ”j 


1 j-1 


a- .esinfiv . ) <$, dx 
2j J k 


( 30 ) 
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where k=l, 2, 3, ... N 


N 

V 

L 


j=i i=i 


n 


rt 2 

( -b. .V . - a. .n m. u . ) u. dx 

lj xi. lj ii. k. 

o JJ i J 1 



(-b. .V . 
lj X J 


a. .Q mu . ) u. dx 
lj j k 



“ N 



+ 

l J 

l (b l3 Vx j i ) - V x j 

' U k 


L j-1 

- i=l 

- 


r ? 

f ' 



| m.a 1 .u. u, dx 

J i lj j. k. 

J 0 J J 1 1 

■ 

l j=i 

l 

- 


+ 


2 

U) 


f R 

ma, . u .u, dx 

i lj J k 


( 31 ) 


where k = 1, 2, 3, ...N 



22 


N n b .V . 

r r r l 1J X J,- 

) ’ ) (a. .u. - ) V , dx 

1 L K lj J i AE. xk i 

j-1 ' i=l 1 


+ [ I! U lj“j - ^ > ***■ ’ 0 


N n 


b. .M . 
4j xj. 


! n. ( vj, ■ -sr 1 > v** 


(a,.$. - 


b. .M . 

-Aljsi ) M dx 

GJ ’ xk QX 


r r l C n 

(a-.v'.' + 77^7- b,.M . 

LJ 2: J EI EI 3j yj 


b 2i M zi } M zk dx f + 

EIj^EIj 3 J i i 


r R C, _ 

! (a_.v'.' + ■ ■■■ b_.M .) 

t 2j j EI x EI 2 3j yj 


UT :: b- .m .) m , } =o 

EI 1 EI 2 2j zj zk 
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j-1 


fR / " 

a 3j W j 


'22 


b„ .M 


EI 1 EI 1 3j yj 


'12 


bn .M_, )M. . dx 


EI^Ij 2j zj yk 


= 0 


(35) 
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5. SELECTION OF FUNCTIONS 

5.1 Basic Functions 

The selection of functions for multiple-load-path blades with fixed- 
conditions for load path members (inboard-ends) is presented in this 
section. The other inboard-end conditions can be incorporated easily 
in a straightforward manner. The CAMRAB program uses nonrotating 
uniform modes for the Galerkin's method; therefore, the same modes are 
employed for the multiple-load-path blades for uniformity. The non- 
rotating modes for uniform cantilever beams are given by: 

1. Bending (Both Flapwise and Chordwise) 
a) Mode shapes 

a f (x) = e r (sina r | - sinh a r |) + 

g (cos a ~ - cosh a §) (36) 

°r r R r R 

where 

e ■ -(sin a - sinh a )/2sin a sinh a 
r r r r r 

g = (cos a + cosh a )/2sin a sinh a 
r r r r r 

r = rth mode shape 

= 1.875, <*2 = 4.694, a g = 7.855, = 10.996 

a 5 = 14.137, a g = 17.279, a ? = 20.420, a g = 23.562 
a 9 = 26.704, a 1Q = 29.845 
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b) Frequencies 



2) Stretching and Twisting 
a) Mode shapes 

a (x) = sin( 6 §) 

r r K 

where 6 = (2r-l) ^ 


(37) 


b) Frequencies 

oo r = 5^ JeA/toR 2 , Stretching 

oo r = 6 r v/GJ/mR 2 , Twisting 


The idealization for a two-load-path blade is shown in Fig. 6. 
The blade is divided into three regions 

Region I: 0 < x < l^, all load path members 

Region II: l <_ x <_ R, blade 

Region III: < x < H, load path members 


For regions I and II, the basic mode functions given by Eqs. (36) 
and (37) are used. It is therefore necessary to develop the clevis 
admissible comparison functions. At this stage the following two options 


are available. 
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1. Satisfy gometric compatibility only and account for the 
force conditions in the Galerkins 1 procedure. 

2. Satisfy both geometric and force conditions. 

Option 1 is selected here because the clevis compatibility functions do not 
depend on the eigenvalues while the clevis admissible comparison functions 
depend on the eigenvalues. 


5.2 Clevis Compatibility Functions 
1. Flapwise Bending 

Assume the functions to be of the following form: 

1,1 ,12,13 

w = a n + a. x + a- x + a_ x (.38; 

r. Or. lr. 2r 3r. 

1 1 1 u 1 

The coefficients in the equation are determined by satisfying the 
displacement and slope continuity at x = 2.^ (see Fig. 6) and 
geometric compatibility at the clevis [Eqs. (13) and (14)]. 

The resulting equations are given by: 


aj; + aj a. + a^ 2,? + a^ t" = w (£..) 
Or. lr , 1 2r. 1 3r . 1 r. 1 

ii i ii 


a. + 2a» 2.. + 3a.. i. = w (2,.) 

lr. 2r. 1 3r. 1 r. 1 

ii ii 


a 0r + V 1 + a 2r ^ + a 3r.‘ 3 * “r (l) + h y.*r U) 

111 1 1 


a} + 2al l + 3a\ l 2 = w'U) 
lr. 2r. 3r. r 

ill 


Arranging the above equations into a matrix form yields: 
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Arranging the above equations into a matrix form yields: 


/» 


— 



- 

r 

s 

1 

a 0r. 

i 


1 

*1 



-l 

w (2- 1 ) 
r . 1 

i 

1 

1 


0 

1 

2*1 

U l 

■ 

w' (£,) 
r . 1 

l 

1 

a o 

2r i 


1 

1 


I 3 


w (£) + h <b r ( 2 -) 
r y i 

3r i , 


0 

1 

21 

\l 2 

J 


w^U) 


where 

subscript r = mode index 

subscript i — load path number index 

superscript 1 = flapwise bending 
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2. Chordwise Bending 

Assume the functions to be of the following form: 


2 2 2 2,2 3 

v = a n + a, x + a- x + a_ x 

r. Or. lr. 2r. 3r. 

1111 l 


(40) 


By employing a similar procedure used for the flapwise bending the 
following equation can be obtained for the chordwise bending clevis 
compatibility function. 







- 

r-i 

r \ 

2 

3 0r. 

l 


1 

*1 


l l 


v<v 

i 

2 

lr . 
i 


0 

1 

2l l 

U l 

. 

> 

V<V 

1 

> 

2 

a o 

2r . 

l 


1 

l 

i 2 



v r (*-) ' h z 

i 

2 

l X- J 


_ 0 

1 

21 

3* 2 . 


v r(*) 


> (41) 


where superscript 2 indicates the chordwise bending. 


3. Twisting 

Assume the function to be of the following form : 


♦ 




l 


(42) 


By satisfying the continuity at x = 2.^ and compatibility at the 
clevis the following relation can be obtained for the coefficients in 
Eq. (42). 
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h 1 1 
b 0r. 

1 

« 

» — 

l 

*1 

-1 

* V 

4> r (<H 1 ) 
i 

i ' 


l 

i 


6(0 

l r i 


The superscript 1 indicates twisting. 


(43) 


4. Stretching 

Assume the functions to be of the following form: 





(44) 


The coefficients in the above equations are computed by 


f h 2 ] 

V 

1 


1 

*1 

-1 

v<h> ) 


’ 3 



■ 


b 2 


1 

l 


u(0 - h w'U) - h v'(i) 

t r i 


- 

- 


z i y i 


(45) 


The superscript 2 in the above equation indicates stretching. 

The treatment for the selection of function is completed so far 
for the deflections, and it is necessary to establish the criteria for 
the forces. For simplicity the functions used for the forces are as 


follows. 
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M = 0 1 

x 


M = w" 

y 


M = v" 
z 


V 

x 


(46) 


For all regions (I, II and III), the basic functions defined in 
Eqs. (36) and (37) are employed while calculating the force functions 
given by Eq. (46). If the following equations are satisfied for 
the force function, then the modified Galerkin' s procedure used in 
CAMRAD would reduce to the standard Galerkin form. 


M = GJd> ' 
x 

M * (EI^cos^P + El^sin^^w" + (E^-EI^ )cosgsin(3v" 
M = (EI,sin 2 3 + EI o cos 2 0)v" + (EI--EI. )cosgsin£w" 

Z I L /I 

V = EAu' 
x 


(47) 


In this situation, it is desirable to use fifth degree polynomials for 
bending degrees-of-freedom and quadratic functions for the twisting 
and stretching degrees-of-freedom for generating the clevis 
compatibility functions. However, the standard Galerkin form is not 
employed for parallel development with the CAMRAD approaches. 
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6 . TENSION CALCULATIONS 

The rotating blade mode shapes calculations are based on the linear 

formulation, and therefore the tension coefficient T appearing in Eqs. 

(2) and (3) should be known in advance; otherwise these equations 

will be nonlinear. Equation (5) can be employed for calculation of 

tension in the blade. The tensions in the load paths are calculated 

by assuming that the load paths are coincident with the blade at the 

clevis, i. .e, h * h =0. The tensions corresponding to this case are 
^i Z i 

given by (Ref. 9): 


fl 


T.(x) 


Q m.xdx + T. (i) 
i x 


(48) 


where 


T.(fc) 


L 


a . 
J 


(49) 


k=l 


n,k 

( n a. ) 
j-1 J 


n,i 

n 

j-i 


a . = a.a. . 
J 12 


. a /a . 
n i 


(50) 


A 


a . = 
i 


dx 


EA. 
o 1 


( 51 ) 
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It is to be noted that the tensions given by Eq. (48) correspond to 


the coincident load path case (h = h =0), and for the noncomcident 

y . z . 

i i 

case (h and/or h ^ 0) the modified Galerkin's method is applied 
y. z. 

J l l 

starting with the coincident load path case. 


Modified Galerkin's Method for Tensions 

The free-body diagram for tension calculations in the load path 
with blade detached is shown in Fig. 7. The equilibrium conditions 
forces at the clevis are given by: 


N 

y v 

L x. 

u i 

i=l 


/■R 


= T(l) = 


2 

Q mxdx 


II 

y v 

L y t 


* 0 


i=l 


II 

)(M + h V - h V ) = 0 

u l ii l i 


i=l 


u 

) (M + h V ) = 0 

L y. z. x. 

i ii 

i=l 


a 

) (M - h V ) = 0 
/ z . y. x. 

u i ^li 


i=l 
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A convenient set of functions are selected for the load path 
members with the following properties. 

1. Satisfy the root boundary conditions of the load 
paths . 

2. Satisfy the clevis compatibility conditions. 

The deflections and forces in the load path members are expanded as 
finite series in terms of the selected functions as shown below 


u . = 

l 


v . = 
l 


w . = 
1 


N 

I 

n=l 

N 

l 

n=l 

N 

l 

n=l 


a. u 
In n . 
x 


a« v 
2n n . 

l 


a., w 
3n n . 

l 


N 

V 


L ^n'V 


n=l 

N 


- I 


b i V 
In xn. 

l 


n=l 

N 


M 


■ l 


b- M 
2n zn . 


n=l 


( 52 ) 
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N 

I b 3 M 

L 3n y° 

n=l 



N 

y b 4 m 

/ 4n xn . 
u 1 

n=l 


The adaptation of Galerkin's method yields the following equations: 


l 

j-1 


1 1 

IJ 


A 

[ b_ .m" . - a (T. w . )' 

3i yj . 3j ij. 
o J J i 


i=l 


+ a, . (Q m. e . xcosg . $ . ) ]w, dx 
4j ii i j . k . 

J J i i 




N 

n 


+ 

l ■ 

) (b_.M . - h b 1 ,V . )w' 

L 3 j yj t z i lj xj. k 

y 


L j=l 

( i=l 

- 


( 53 ) 
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N f n 

mi 


n=l 


[b ? .M . 
o 2j ZJ i 


a-,.(T.v!)' + a, . (Q m . e . xsinB . $ . ) 
2j i j. 4j ii p i y j 


i=l 


2 2 
+ a. .Q me . e . sing . 0 . - a».Q m.v. ]v. dx 

4 j i i i Y j. 2 j l k. 


N n 

) i ) (-b_ .M 1 . + a . T.v! 

L \ L 2j Zj 2j l J, 

_ j=l ' i=l 


- a. . n m. e . xsin$ .<fi. )v, > 
4j ii i k 


4 x=e 


N 

f ? 


l \ 

1 (b 2 j’ M zV h y. b l j V x j . )v k 

> 

_ j =1 

i=l 



= o 


4 x=£ 


(54) 


N 

V 

L 

j=i 


c 

n 


i=i 


-b. . m' . + e ,x( -a- . cos3 . v! 

4j xj, ii 2j i j 


+ a_.sin3.wl ) + a_.fl m.e.sinp.v. 
3j i 2j i i p i j. 


+ a. .Q m.(k -k )cos 23 .<f. 
4j l m_ m. ii. 

2.1. J l 

i l 


r <)> k dx 

i 
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N n 


) ■ ) [b. .M . + h ( -b- .M 1 . + a-.T.w'. 

L L 4 J xj. y i 3j yj. 3j i j. 


j=l v i-1 


+ a, .Q m . e . xcosjS . <f . ) - h (-b_.M'. 
4 J i i i Jj_ 2j zj. 


, 2 

+ a_.T.v. - a . .Q m. e . xsinfi . 0 . )] 0, r > = 0 

2j l j. 4j l l H i"j. T k 

x=£ 


N f n 


y y (-b. .V 1 . - ai.fl^m.u, )u. dx 

L ) L J Q ij xj. ij i j,. k i 



b. -V . u, 
Ij xj. k 


« T(£)u k (l) 


N n 


1 I 


b lj V xj. 

-EA-^ )V xk. dx - 0 
1 1 


N n 


b. .M . 

4 J xj. 


(a,,*. — ) M , dx } = 0 


/ V 4 J ? J. GJ. ' xk. 
J o 1 1 


N n 


T ■ f ^ »v» ■ 

L L J o 2] J. EI 1 EI 2 3j yj. 


n T ‘ nj b. .M . )M , dx f = 0 
EI 1 EI 2 2 j zj. zk. 
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N n 

l l 

j=l ^ i=l 


r ( - _i2L 

J o a 3j W j i - EI lE I 2 k '3j“yj 1 


b„ .M 


'12 


b_ ,M . )M , dx 


2j zj i yk i 


= 0 


(60) 


where c^, c^, c 22 are defined in Eq. (35). The solution of 

Eqs. (53) to (60) yields T^(2.) [which is (£.)], and from Eq. (48) 

i 

the tension coefficients in the load path members can be computed. At 
the most, two or three iterations may be needed starting with the coinci- 
dent load path tensions to solve the nonlinear problem. The functions 
developed in Chapter 5 are used to solve the present problem a-lso. 
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7 ■ SOFTWARE DEVELOPMENT 


7.1 Equations of Motion 

The basic equations of motion given by Eqs. (18) to (25) are 
nondimensionalized as shown below. 


V 

du x 
dx EA 


(61) 


where u = u/R; x = x/R; EA = EAR 2 /EI t 


V = V R /EI n 
xx R 


,r M 
d<fr x 

dx GJ 


(62) 


where 


GJ = G J/EI d ; M = M R/EI d 

K X X K 


12 


.2- -C. 

d v 

.-2 ~ EI 1 EI_ 
dx 12 


M + M 

y EI^EI^ z 


(63) 


where = (E^-EI^cosjjsinji 


C n = El^cos^fj + El^sin^p 


EIj. = EI 1 /EI r ; EI 2 = EI 2 /EI r ; v = v/R 


\ = V /E V = M 3 R/EI R 
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where 


where 


where 


where 


-ill- 5 

dx £ I 1 EI 2 y ei 1 ei 2 2 

— 2 - — 2 - 
C 22 = EI i sin 13 + E l2 cos & 

w = w/R 


( 64 ) 


-dv 
x 

dx 


=. 2 — 2 — 
- S2 mu = uj mu 


( 65 ) 


— 9 9 A 

Q =12 m^R /EI^; m = m/m^ 
uj 2 = u^m^R^/EI^ 

d ^x -2 -2 -2 

~ + 12 m(k -k )cos 234 > 

dx m 2 ra l 

“ 2 ”"” _ _ « 2 -— - - 

- 12 mex(-singv' cosgw') + 12 mesinpv 




- 2 -- 2 - - 2 -- -- - 2 - - 
uj ink di - uj mesingv + uj mecospw 


k 2 /R 2 ; k = k 2 /R 2 ; 
m 2 m 2 


‘ iX * ;-«/*, ( )' -£ 


- (T ^) 1 - (Q 2 mexcos|3$) 1 
dx 

- 2 -- , - 2 -- -- 
= uj mw + uj mecos(3$ 



mxdx 


( 66 ) 


( 67 ) 
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— -rr - (T ~) ' + (^mexsin^)' 
dx dx 

- 2 -- -- - 2 -- 
+ Q mesinpd - Q mv 

- 2 -- - 2 -- 

= (jj mv - oj mesinfSd (68) 

Nonuniform properties are allowed in the computer program and maximum 
number of stations allowed for inputting the properties are 11 for each 
load path member and 51 for the blade. The linear variation of 
properties is assumed for the intermediate locations. 

7.2 Assumed Functions 

The selected basic functions are given by Eqs. (36) and (37). 

These functions are nondimensionalized and differentiated with respect 
to x as shown below. 

1. Bending (Both Flapwise and Chordwise) 
a) Deflections and Derivatives: 

a (x) = e (sin a x - sinh a x) + 
r r r r 

g (cos a x - cosh a x) (69) 

®r r r 

a'(x) = a e (cos a x - cosh a x) + 
l r r r r 

a g (-sin a x - sinh a x) 
rr r r 


( 70 ) 
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b) Moments and Derivatives: 


M (x) = a^e (-sin ax- sinh a x) + 
r r r r r 


a^g (-cos a x - cosh a x) 
r°r r r 


(71) 


M (x) = a^e (-cos ax- cosh a x) + 
r r r r r 


a^g (sin a x - sinh a x) 
rr r r 


(72) 


where 


a = a , g = g , e = e 
r r r r r r 


and are defined in Eq. (36). . 


2. Twisting/Stretching 
a) Deflections 


a^_(x) = sin(6 r x) 


(73) 


b) Moment/Force 


M (x) = 6 cos(6 x) 
r r r 


(74) 


where 6^ = (2r-l) ^ 

The clevis compatibility functions are completed assuming 
= 0.9 and the coefficients for the bending deflection functions 
are stored in the triple subscripted variable A(I,J,K) with subroutine 
titled GFCTS. The notation for the arguments is as follows: 



I = Number of selected functions 
J = i - a Q 

J = 2 -*• 

J = 3 - a 2 

J = 4 •> a 3 

where a^, a^, are coefficients of the cubic polynomial 

K = 1 -*■ Flapwise bending, Load path 1 

K = 2 -*■ Chordwise bending, Load path 1 

K = 3 -+ Flapwise bending. Load path 2 

K = 4 ■> Chordwise bending. Load path 2 


The coefficients for the twisting/stretching deflection functions 
are stored in the triple subscripted variable B(I,J,K) in the 
subroutine GFCTS . The notation for the arguments is as follows: 


where 


I = 
J = 
J = 

V 

K = 
K = 
K = 
K = 


Number of selected function 


b^ are coefficients of the linear function 

1 -* Twisting, Load path 1 

2 -*■ Stretching, Load path 1 

3 -*■ Twisting, Load path 2 

4 -> Stretching, Load path 2 
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7.3 Galerkin's Equations 

The Galerkin's equations to solve the eigenvalues problem are given 
by Eqs . (28) to (31). To avoid the numerical differentiation of several 
products of the nonuniform blade properties, the above equations are 
rewritten as shown below. 

(b_ M - a_.D113.w! + a..D12.$! )w' dk f 

3 j yj t 3 J i J t i J. k i 
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b- .M . ]w. > 

3 j yj k 


x= l 
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+ u 2 


l J 
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(D21a_.w. + a. ,D23<J> . )w. dx 
3 J J J k 


(75) 
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*1 
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o J J i J J 11 


i=l 


46 


N 

r 

f 

< 


L 

. 

=1 

t 

N 

r 

r 

/• 



L 

. 

=i 

t 


t <b 2j M Ij ' a 2j D213v j + a 4 j D214, *j> v k dx 


( a , . Q 2 D2S<(. . - a, ,0 2 D21.v.)v, dx 
J 2j 1 j k 


N 

" 

" 

- 

+ < l 


L 'V.j, • VuV 

• - b 9 .M . 
2j zj 

t—i 

II 

**-) 

- 

1 i-1 

- 


V k > 


x=£ 


N 


= 0) 


2 


HU ( D 11 i a 2 jV - a 4 j D 15 i^j 


l l 


j-1 ' i-1 


2 V ! 

+ U) L}) 

j=l 


rR 

(D21a-.v. - a . .D250 . )v, dx }■ 
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(76) 


I ' y ( ~b ,M 1 . -Dl2 a 9 .vl + D14 a w! 

L L J 0 4j xj. i 2j J± i 3] j. 

J-1 i-1 


-2 

+ a- .Q D15.v. + a. .D16.0. )$, dx 
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+ < 
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Dil3 » T 

i = 1 -> Load Path Properties 
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2 Blade Properties 
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7.4 Eigenvalue Problem 

Eqs. (75) to (82) can be written into a set of matrix equations 
of the following form for k = 1 to N. 


[A 1 ] 

NX4N 


3 

4NX1 


w [b l j 

N.X2N 


4 

2NX1 


(83) 



r > 
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l 4 J 

2NX1 


(84) 


[a 3 ] 

NX6N 
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- w 2 [B 3 J 
NX3N 



3NX1 


(85) 
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CA 4 1 

NX2N 



“ 2 Ib 4 J 

NXN 



( 86 ) 



(87) 



( 88 ) 


[A ? ] 

NX3N 



(89) 


[a 8 ] 

NX3N 



(90) 


The above set of equations correspond to Eqs. (75) to (82) respectively. 
Eq. (83) can be partitioned into the following form: 
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[ A u] | a 12 j 


- - 3 - t + 


l ^ i 


U13] , [A u ] 


(91) 


= . IB U ] | [B 12 ] 


By employing a similar partition scheme for Eqs. (84) to (86), the 
four equations (83 to 86) can be combined into the following form 
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[E 1 3 = 
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B 11 
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B,-, 

0 

0 

0 


Eqs. (87) to (90) can be combined into a matrix equation of the 
following form using a similar partition scheme employed in Eq. 
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Substituting Eq. 

(3) 

into 

Eq. (92) yields 
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81 


(91). 

(93) 
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(tDj - u) Z [I] ) 



( 94 ) 


where 


[D] = [E 3 ] _1 ( fE 1 ] + [E 2 J [E 4 D 



56 


REFERENCES 


1. Fenaughty, R. R. and Noehren, W. L., "Composite Bearingless Tail 
Rotor for UTTAS," Journal of the American Helicopter Society, 
vol. 22, (3), July 1977. 

2. Staley, J. A., Gabel, R. and MacDonald, H. I., "Full Scale Ground 
and Air Resonance Testing of the Army-Boeing Vertol Bearingless 
Main Rotor," Paper No. 79-23, American Helicopter Society 

35th Annual National Forum, Washington, D.C., May 1979. 

3. Dixon, P. G. C. and Bishop, H. E. , "The Bearingless Main Rotor," 
Journal of the American Helicopter Society, vol. 25 (3), July 1980. 

4. Bousman, W. G., Ormiston, R. A. and Mirick, P, H., "Design 
Considerations for Bearingless Main Rotor Hubs," American Helicopter 
Society 39th Annual National Forum, St. Louis, MO, May 1983. 

5. Johnson, Wayne, "A Comprehensive Analytical Model of Rotorcraft 
Aerodynamics and Dynamics," NASA TM 81182, June 1980. 

6. Lang, K. W. and Nemat Nasser, S., "An Approach for Estimating 
Vibration Characteristics of Nonuniform Rotor Blades," AIAA Journal, 
vol. 17, No. 9, Sep. 1979. 

7. Houbolt, J. C. and Brooks, G. W., "Differential Equations of Motion 
for Combined Flapwise Bending, Chordwise Bending, and Torsion of 
Twisted Nonuniform Rotor Blades," NACA Report 1346, 1958. 

8. Dixon, P. G. C., "Design, Development and Flight Demonstration of 
the Loads and Stability Characteristics of Bearingless Main Rotor," 
USA AVRADCOM-TR- 80 -D3 , 1980. 

9. Murthy, V. R., and Joshi, A. M., "Free Vibration Characteristics 
of Multiple-Load-Path Blades by the Transfer Matrix Method," 

Journal of the American Helicopter Society, vol. 31, (4), 

pp. 43-50, Oct. 1986. 



PROGRAM MOGALK 


C ***************************************************************** 

C DYNAMIC ANALYSIS OF MULTIPLE -LOAD -PATH BEARINGLESS ROTOR BLADES 

C MODIFIED GALERKIN' S METHOD COUPLED FLAPWISE BENDING, CHORDWISE 

C BENDING, TWISTING, AND STRECHING MOTION OF PRETWISTED NONUNIFORM 

C ROTATING BLADES . 

Q ★★★★★**★★*******★*;*:★★**★*★★*★;*★★****★******************★★★*★***** 

REAL MASS , MASS1 (2,11) ,MASS2 (51) ,KM1S1 (2,11) ,KM1S2 (51) ,KM2S2 (51) , 
1KM2S1 (2,11) 

DIMENSION STA1 (11) , STA2 (51) , sta (51) ,EIY1 (2,11) ,EIY2 (51) ,EIZ1 (2, 
111) ,EIZ2 (51) , EA1 (2,11) , EA2 (51) , BETA1 (2,11) , BETA2 (51) , GJ1 (2,11), 
2GJ2 (51) ,E1 (2, 11) ,E2 (51) ,FD (2,2) ,D21 (51) ,D22 (51) ,D23 (51) ,D24 (51) , 
3D25 (51) ,D2 6 (51) ,D27 (51) ,D28 (51) ,D29 (51) ,D210 (51) ,D211 (51) , 

4D212 (51) ,D213 (51) ,D11 (2, 11) ,D12 (2,11),D13(2,11) ,D14 (2,11) , 

5D15 (2,11) ,D1 6 (2,11) ,D17 (2,11) ,D18 (2, 11) , D19 (2, 11) ,D110 (2,11), 
6D111 (2, 11) ,D112 (2,11) ,D113 (2, 11) ,BB(3) ,AT (2) ,DS (20,20) , 

7A (5,4,4) ,B(5,2,4) ,D214 (51) ,FRE(51) 

COMPLEX DSC (20, 20) , LAMDA (20) ,VECT(20,20) ,HL (20, 20), HI (20, 20), 

1CNT (20) , COLM (20) 

INTEGER INTH (20,2) 

COMMON /XM1 /NPATH 

COMMON /XM2 / SPAN1 , SPAN2 , SPAN3 , FD , OMEGAN 

COMMON/XM3/D11 , D12,D13,D14,D15,D16,D17,D18,D19,D110,D111,D112, 
1D113,D21,D22,D23,D24,D25,D26,D27,D28,D29,D210,D211,D212,D213 
COMMON/XM4/STA, NS 
COMMON / XM5 / HS 1 , H S 2 
COMMON/XGFCT1/A, B 
COMMON/XGPR1/DS 
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c 

C THIS SECTION READS THE DATA OF THE SYSTEM 

C 

C 

READ (20, 100) NPATH, NS1,NS2, MODES, ITT, NV 
READ (20,105) SPAN1, SPAN2, OMEGA 
READ (20,105) (STA1 (I) , I=1,NS1) 

READ (20, 105) (STA2 (I) , 1=1, NS2) 

DO 400 1=1, NPATH 

READ (20,105) (MASS1 (I , J) , J=1 , NS1 ) 

READ (20,105) (EIY1 (I, J) , J=1,NS1) 

READ (20,105) (EIZ1 (I, J) ,J=1,NS1) 

READ (20,105) (GJ1 ( I , J) , J=1 , NS1 ) 

READ (20,105) (El (I , J) , J=1 , NS1 ) 

READ (20, 105) (BETA1 ( I , J) , J=1 , NS1 ) 

READ (20,105) (KM1S1 (I, J) , J=l, NS1) 

READ (20, 105) (KM2S1 ( I , J) , J=1 , NS1 ) 

READ (20,105) (EA1 (I, J) , J=1,NS1) 

400 CONTINUE 

READ (20,105) (MASS2 (I) , 1=1, NS2) 

READ (20,105) (EIY2 (I) , 1=1, NS2) 

READ (20,105) (EIZ2 (I) ,I=1,NS2) 

READ (20,105) (GJ2 (I) , 1=1, NS2) 

READ (20,105) (E2 (I) , 1=1, NS2) 

READ (20, 105) (BETA2 (I) , 1=1, NS2) 

READ (20,105) (KM1S2 ( I ) , 1=1 , NS2 ) 

READ (20,105) (KM2S2 (I ) , 1=1 , NS2 ) 

READ (20,105) (EA2 (I) , 1=1, NS2) 
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DO 401 1=1 , NPATH 
DO 401 J-1,2 
401 FD ( I , J) =0 . 0 

IF (NPATH. GT.l) READ (20, 105) ( (FD (I , J) , J=1 , 2 ) , 1=1 , NPATH) 


THIS SECTION 

PRINTS THE DATA OF THE SYSTEM 

SPAN=SPAN1 * SPAN2 
WRITE (22,200) 
WRITE (22,205) 

WRITE (22,347) 

MODES 

WRITE (22,371) 

NV 

WRITE (22,225) 

SPAN 

WRITE (22,230) 

SPAN1 

WRITE (22,235) 

SPAN2 

WRITE (22,240) 

OMEGA 

WRITE (22,245) 

NS1 

WRITE (22,250) 

NS2 

WRITE (22,315) 

NPATH 

WRITE (22, 316) 

ITT 


IF (OMEGA. LE.O) ITT=1 
DO 405 1=1, NPATH 

IF (I.EQ.l) WRITE (22, 255) 

IF (I.EQ.2) WRITE (22, 260) 

WRITE (22,270) 

WRITE (22,275) (STA1 ( J) , J=l, NS1 ) 
WRITE (22,280) 
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WRITE (22,275) (MAS SI (I, J) , J=1,NS1) 

WRITE (22,285) 

WRITE (22,275) (EIY1 (I, J) , J=1,NS1) 

WRITE (22,350) 

WRITE (22,275) (EIZ1 (I, J) , J=1,NS1) 

WRITE (22,352) 

WRITE (22,275) (GJ1 (I, J) , J=1,NS1) 

WRITE (22,354) 

WRITE (22, 275) (El (I , J) , J=1 , NS1) 

WRITE (22,356) 

WRITE (22,275) (BETA1 (I, J) , J=1,NS1) 

WRITE (22,358) 

WRITE (22,275) (KM1S1 (I, J) , J=1,NS1) 

WRITE (22,360) 

WRITE (22,275) (KM2S1 (I, J) , J=1,NS1) 

WRITE (22, 286) 

WRITE (22, 275) (EA1 (I, J) , J-1,NS1) 

IF (NPATH.GT.l) WRITE (22 , 321 ) I , FD (I , 1 ) , FD (I , 2 ) 
405 CONTINUE 

WRITE (22,290) 

WRITE (22,270) 

WRITE (22,275) (STA2 (I) ,I=1,NS2) 

WRITE (22,280) 

WRITE (22,275) (MASS2 (I) , 1-1, NS2) 

WRITE (22,285) 

WRITE (22,275) (EIY2 (I) , 1-1 , NS2) 

WRITE (22, 350) 

WRITE (22,275) (EIZ2 (I) ,I=1,NS2) 
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WRITE (22,352) 

WRITE (22,275) (GJ2 (I) , 1=1, NS2) 
WRITE (22,354) 

WRITE (22,275) (E2 (I) , 1=1, NS2) 

WRITE (22,356) 

WRITE (22,275) (BETA2 (I ) , 1=1 , NS2) 
WRITE (22,358) 

WRITE (22,275) (KM1S2 (I ) , 1=1 , NS2 ) 
WRITE (22,360) 

WRITE (22, 275) (KM2S2 (I) , 1=1, NS 2) 
WRITE (22,286) 

WRITE (22,275) (EA2 (I) , 1=1, NS2) 
NS=NS1 

DO 406 J=1 , NS 
STA(J) =STA1 (J) 

406 CONTINUE 

HH=SPAN1/10 . 0 
DO 418 1=1 , NPATH 
DO 410 J=1 , NS1 

D21 (J)=MASS1 (I, J) 

D22 ( J) =EIY1 (I, J) 

D23 ( J) =EIZ1 (I, J) 

D24 ( J) =GJ1 (I, J) 

D25 ( J) =E1 (I, J) 

D26 (J) =BETA1 (I, J) 

D27 (J)=KM1S1 (I, J) 

D28 (J) =KM2S1 (I, J) 

D29 (J) =EA1 (I, J) 



410 CONTINUE 


CALL INTPOL (11, D21,HH) 

CALL INTPOL (11, D22, HH) 

CALL INTPOL (11, D23,HH) 

CALL INTPOL ( 11 , D24 , HH) 

CALL INTPOL (11, D25,HH) 

CALL INTPOL (11, D2 6, HH) 

CALL INTPOL (11, D27,HH) 

CALL INTPOL (11, D2 8, HH) 

CALL INTPOL ( 11, D2 9, HH) 

DO 415 J-1,11 

MASS1 (I, J)=D21 (J) 

EIY1 (I, J)«D22 (J) 

EIZ1 (I, J)-D23 (J) 

GJ1 (I, J) =D24 (J) 

El (I, J)-D25 (J) 

BETA1 (I, J) =D26 ( J) 

KM1S1 (I, JJ-D27 (J) 

KM2S1 (I, J)«D28 (J) 

EA1 (I, J) =D29 ( J) 

415 CONTINUE 

IF (I.EQ.l) WRITE (22, 255) 

IF (I.EQ.2) WRITE(22, 260) 
WRITE (22,280) 

WRITE (22,295) 

WRITE (22,275) (D21 ( J) , J-l, 11) 

WRITE (22,285) 

WRITE (22,275) (D22 ( J) , J-l, 11) 
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WRITE (22,350) 

WRITE(22,275) (D23 (J) , J-l, 11) 

WRITE (22,352) 

WRITE (22,275) (D24 ( J) , J-l, 11) 

WRITE (22,354) 

WRITE (22,275) (D25 (J) , J-l, 11) 

WRITE (22,356) 

WRITE (22,275) (D26 { J) , J-l, 11) 

WRITE (22,358) 

WRITE (22,275) (D27 ( J) , J-l, 11) 

WRITE (22,360) 

WRITE (22,275) (D28 ( J) , J-l, 11) 

WRITE (22,286) 

WRITE (22,275) (D29 ( J) , J-l, 11) 

418 CONTINUE 
NS-NS2 

DO 419 J-l, NS 
STA ( J) =STA2 (J) 

419 continue 
HH-SPAN2/100.0 

CALL INTPOL (51 , MASS2 , HH) 

CALL INTPOL (51, EIY2,HH) 

CALL INTPOL (51, EIZ2,HH) 

CALL INTPOL (51, GJ2,HH) 

CALL INTPOL (51, E2,HH) 

CALL INTPOL (51, BETA2,HH) 

CALL INTPOL (51, KM1S2,HH) 

CALL INTPOL (51, KM2S2,HH) 
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CALL INTPOL (51 , EA2 , HH) 

CALL INTPOL (51, CH,HH) 

CALL INTPOL (51 , QE, HH) 

WRITE (22,300) 

WRITE (22, 280) 

WRITE (22,275) (MASS 2 ( J) , J=l, 51) 

WRITE (22,285) 

WRITE (22, 275) (EIY2 ( J) , J=l, 51) 

WRITE (22,350) 

WRITE (22,275) (EIZ2 ( J) , J-l, 51) 

WRITE (22,352) 

WRITE (22, 275) (GJ2 ( J) , J-l , 51 ) 

WRITE (22,354) 

WRITE (22,275) (E2 ( J) , J-l, 51) 

WRITE (22,356) 

WRITE (22, 275) (BETA2 (J) , J-l, 51) 

WRITE (22,358) 

WRITE (22, 275) (KM1S2 ( J) , J-l, 51) 

WRITE (22,360) 

WRITE (22, 275) (KM2S2 ( J) , J-l , 51) 

WRITE (22,286) 

WRITE (22, 275) (EA2 ( J) , J-l , 51 ) 

WRITE (22,205) 

C 

C THIS SECTION NOND IMEN S IONALI ZED THE DATA 

C 

C 


PI-4. 0*ATAN (1.0) 
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OMEGA=OMEGA*P I / 3 0 . 0 
OME GAS=OME GA* OME GA 
SPANS=SPAN* SPAN 
EIY=EIY2 (1) 

MASS=MASS2 (1) 

FACT=SQRT (SPANS*SPANS*MASS/EIY) 

OME GAN =OME GAS * FACT* FACT 
FA=SPANS/EIY 
SPAN1=SPAN1 /SPAN 
SPAN2=SPAN2 /SPAN 
HS1=SPAN1/10 . 0 
HS2=SPAN2 / 5 0 . 0 
SPAN3=0.7*SPAN1 
DO 420 1=1 , NPATH 
DO 420 J-1,11 

BETA1 (I, J)=BETA1 (I, J) *PI/180.0 

KM1S1 (I, J)=KM1S1 (I, J) / (MASS1 (I, J) *SPANS) 

KM2S1 (I, J)-(KM2S1 (I, J) /MASS1 (I, J) +E1 (I, J) *E1 (I , J) ) /SPANS 
El (I, J)=E1 (I, J) /SPAN 
MAS SI (I, J)=MASS1 (I, J) /SPAN 
EIY1 (I, J)-EIY1 (I, J) /EIY 
EIZ1 (I, J) =EIZ1 (I, J) /EIY 
GJ1 (I, J)=GJ1 ( I , J) /EIY 
EA1 (I, J)=EA1 (I, J) *FA 
420 CONTINUE 


DO 4201 J=1 ,51 

BETA2 ( J) = (BETA2 (J) ) *PI/180.0 
KM1S2 ( J) =KM1S2 ( J) / (MASS2 ( J) *SPANS) 
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KM2S2 ( J) = (KM2S2 (J) /MASS2 (J) +E2 (J) *E2 (J) ) /SPANS 

E2 (J)=E2 (J) /SPAN 

EIY2 (J)=EIY2 (J) /EIY 

EIZ2 ( J) =EIZ2 (J) /EIY 

GJ2 ( J) =GJ2 (J) /EIY 

EA2 ( J) =FA*EA2 (J) 

4201 CONTINUE 

C 

C THIS SECTION CALCULATES THE COEFFICIENTS IN THE GALERKIN' S EQUATIONS 

C 

c 

x=o . 0 

DO 430 1=1 , NPATH 
DO 430 J-1,11 
TC1=0MEGA*MASS1 (I, J) 

C=COS (BETA1 (I, J) ) 

S=SIN (BETA1 (I, J) ) 

Dll (I, J)=MASS1 (I, J) 

D12 (I, J) =TC1 *E1 (I, J) *X 

D14 (I, J)=D1 2 (I, J) *S 

D12 (I, J)=D12 (I, J) *C 

D13 (I, J)=MASS1 (I, J) *E1 (I, J) 

D15 (I/J)=D13 (I, J) *S 
D13 (I, J)=D13 (I, J) *C 

D16 (I, J) =TC1 * (KM2S1 (I, J) -KM1S1 (I, J) ) * (C*C-S*S) 

D17 (I, JJ-MASS1 (I, J) * (KM2S1 (I, JJ+KM1S1 (I, J) ) 

D18(I f J)=l. 0/EA1 (I, J) 

D19 (I, J) =1 . 0/GJ1 (I, J) 



TC2=EIY1 (I, J) *EIZ1 (I, J) 

DUO (I, J) = (EIY1 (I, J) *C*C+EIZ1 (I, J) *S*S) /TC2 
Dill (I , J) = (EIZ1 (1/ J) -EIY1 (I, J) ) *C*S/TC2 
D112 (I, J) = (EIY1 (I, J) *S*S+EIZ1 (I, J) *C*C) /TC2 
X=X+HS1 
430 CONTINUE 
X-SPAN1 
DO 440 J-1,51 
TCl=OMEGAN*MASS2 ( J) 

C=COS (BETA2 ( J) ) 

S=SIN (BETA2 (J) ) 

D21 ( J) =MASS2 (J) 

D22 (J) =TC1*E2 ( J) *X 
D24 (J) =D22 ( J) *S 
D22 ( J) =D22 (J) *C 
D23 ( J) =MASS2 ( J) *E2 ( J) 

D25 (J)=D23 (J) *S 
D23 ( J) =D23 (J) *C 

D26 ( J) =TC1* (KM2S2 ( J) -KM1S2 ( J) ) * (C*C-S*S) 

D27 ( J) =MASS2 ( J) * (KM2S2 ( J) +KM1S2 ( J) ) 

D28 (J) =1 . 0/EA2 ( J) 

D29 ( J) =1 . 0/GJ2 (J) 

TC2=EIY2 ( J) *EIZ2 (J) 

D210 ( J) = (EIY2 (J) *C*C+EIZ2 (J) *S*S) /TC2 
D211 ( J) = (EIZ2 ( J) -EIY2 (J) ) *C*S/TC2 
D212 ( J) — (EIY2 (J) *S*S+EIZ2 (J) *C*C) /TC2 
X=X+HS2 


440 CONTINUE 
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X=SPAN1 
H=HS2/24 . 0 
DO 450 1=1,51 
D214 (I)=MASS2 (I) *H*X 
X=X+HS2 
450 CONTINUE 

D213 (51 ) =0 . 0 

D213 (50) = (9 . 0*D214 (51 ) +19 . 0*D214 (50 ) -5 . 0*D214 (49)+D214 (48) ) 
l*OMEGAN 
DO 455 1=2,49 
J-51-I 

D213 ( JJ-D213 (J+l) + (-D214 (J+2) +13.0* (D214 (J+l)+D214 ( J) ) - 
1D214 ( J-l) ) *OMEGAN 
455 CONTINUE 

D213 (1) =D213 (2) + (D214 (4) -5. 0*D214 (3 ) +1 9 . 0*D214 (2) +9.0*D214 (1) ) 

1 *OMEGAN 

GO TO (460, 465) ,NPATH 
460 BB (1) =D213 (1) 

GO TO 480 
465 H=SPAN1/ 24.0 

DO 475 J=1 , NPATH 

AT ( J) = (9 . 0*D18 ( J, 1) +19 . 0*D18 ( J, 2) -5 . 0*D18 ( J, 3) +D18 (J,4) ) *H 
DO 470 K-1,9 

AT ( J) =AT ( J) +H* (-D18(J,K)+13.0*D18 ( J, K+l ) +13 . 0*D18 (J,K+2) 

1-D18 ( J, K+3 ) ) 

470 CONTINUE 

AT (J) =AT ( J) +H* (D18 (J, 8) -5. 0*D18 ( J, 9) +19 . 0*D18 (J, 10) +9. 0*D18 (J, 11) ) 


475 CONTINUE 
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BB ( 1 ) =AT ( 1 ) *D213 (1) / (AT(1)+AT(2) ) 

BB ( 2 ) =AT ( 2 ) *D213(1) / (AT(1)+AT(2) } 

480 DO 495 JJ=1, NPATH 
X=0 . 0 

DO 485 1=1,11 
D214 (I) -MASS1 ( JJ, I) *H*X 
X=X+HS1 
485 CONTINUE 

D113 (JJ,11)=0.0 

D113 (JJ,10)=(9.0*D214 (11) +19 . 0*D214 (10) -5 . 0*D214 ( 9) +D214 (8) ) 
l*OMEGAN 
DO 490 1=2,9 
J=ll-I 

D113 (JJ, J)=D113 (JJ, J+l) +C-D214 (J+2)+13.0* (D214 (J+l)+D214 (J) 
1-D214 ( J-l) ) *OMEGAN 
490 CONTINUE 

D113 ( JJ, 1) =D113 ( JJ, 2) + (D214 (4)-5.0*D214 (3) +19 . 0*D214 (2) 

1+9. 0*D214 (1) ) *OMEGAN 
495 CONTINUE 

DO 500 1=1 , NPATH 
DO 500 J=l,ll 

500 D113 (I,J)=D113 (I, J) +BB (I) 

IF (NPATH. GT. 1 .AND. ITT. GT . 1) CALL TITER (ITT) 

C 

C THIS SECTION CALCULATES GALERKIN'S FUNCTIONS AND WRITES THE 

C COEFFICIENTS OF THE CLEVIS COMPATIBILITY POLYNOMIALS 

C 


C 
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CALL GFCTS (MODES) 

Nl=2 *NPATH 
WRITE (22,361) 

DO 505 K-l, MODES 
WRITE (22,366) 

DO 505 1=1, Nl, 2 
IF (I.EQ.l) WRITE (22, 362) 

IF (I.EQ.2) WRITE (22, 363) 

WRITE (22,364) 

WRITE (22,2 75) (A(K, J, I) , J-l, 4) 

WRITE (22,365) 

WRITE (22, 275) (A (K, J, 1+1) , J-l, 4) 

WRITE (22,367) 

WRITE (22,275) (B (K, J, I) , J=l,‘2) 

WRITE (22,368) 

WRITE (22,275) (B (K, J, 1+1) , J-l, 2) 

505 CONTINUE 

C 

C THIS SECTION CALCULATES THE EIGENVALUES AND EIGENVECTORS 

C 

C 

N4=4*NPATH 
DO 510 1=1, N4 
DO 510 J— 1 , N4 

510 DSC (I, J) =CMPLX (DS(I,J) ,0.0) 

INTH (1,1) =N4 
INTH (2, 1) -NV 

CALL EECM (DSC, LAMDA, VECT, HL, HI, CNT, COLM, INTH) 
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WRITE (22,200) 

WRITE (22,205) 

WRITE (22,369) 

WRITE (22,205) 

WRITE (22,275) (LAMDA(J) , J=1,N4) 

DO 515 J=1,N4 

515 FRE ( J) =SQRT (REAL (LAMDA( J) ) ) 

WRITE (22,370) 

WRITE (22, 275) (FRE ( J) , J=1 , N4 ) 

IF (NV.EQ.O) GO TO 600 
WRITE (22,371) 

DO 520 1=1, N4 

520 WRITE (22,275) (VECT (I , J) , J-l , NV) 

WRITE (22,372) 

DO 530 J=1,NV 
AMAX=REAL (VECT ( I , J) ) 

DO 525 1=2, N4 

525 IF (ABS (AMAX) . LT . ABS (REAL (VECT ( I , J) ) ) ) AMAX=REAL (VECT ( I , J) ) 
DO 530 1=1, N4 

530 VECT (I, J) =VECT (I, J) /AMAX 
DO 535 I=N4 

535 WRITE (22, 275) (VECT (I , J) , J-l , NV) 

WRITE (22,205) 

C 

C FORMATS 

C 


C 



100 FORMAT (515) 

105 FORMAT (5E14 . 7) 
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200 FORMAT (1H1) 

2 05 FORMAT (//2X/ ' *****************************************^ f 

£ t ★★**********'**************f ) 

225 FORMAT (//5X, 'RADIUS OF THE ROTOR (INCHES) =',E14.7) 

230 FORMAT ( / / 5X, ' LENGTH OFF INBOARD SEGMENTS (IN) -' , E14 . 7) 

235 FORMAT ( //5X, ' LENGTH OF THE BLADE (INCHES) -',E14.7) 

240 FORMAT (//5X, 'ROTATIONAL SPEED (RPM) =',E14.7) 

245 FORMAT (//5X, 'NUMBER OF DATA POINTS FOR INBOARD SEGMENTS 3 ' , 15) 
250 FORMAT (//5X, 'NUMBER OF DATA POINTS FOR BLADE -',15) 

255 FORMAT (//5X, 'PROPERTIES OF THE FIRST LOAD PATH') 

260 FORMAT (//5X, 'PROPERTIES OF THE SECOND LOAD PATH') 

270 FORMAT (//5X, 'DATA POINT LOCATIONS IN INCHES') 

275 FORMAT (4 ( 6X, E14 . 7 ) ) 

280 FORMAT (//5X, 'MASS PER UNIT LENGTH (LB-SEC . **2/IN**2 )' ) 

285 FORMAT (//5X, 'FLAPWISE BENDING STIFFNESS (LB-IN**2)') 

286 FORMAT (//5X, 'AXIAL STIFFNESS (LB) ') 

290 FORMAT (//5X, 'PROPERTIES OF THE BLADE') 

295 FORMAT (//5X, ' INTERPOLATED VALUES FOR INBOARD SEGMENTS, 11', 

$' EQUIDISTANT VALUES') 

300 FORMAT (//5X, ' INTERPOLATED VALUES FOR THE BLADE, 51', 

$' EQUIDISTANT VALUES') 

315 FORMAT (//5X, 'MUNBER OF LOAD PATHS') 

316 FORMAT (//5X, 'NUMBER OF FLNSION ITERATIONS -',15) 

321 FORMAT (//5X, 'Y, Z DISTANCES BETWEEN LOADPATH NO: ',12, 

$/5X, ' AND THE BLADE ( IN) ARE 3 ' , F6 . 2 , 5X, F6 . 2 ) 

347 FORMAT (//5X, 'NUMBER OF MODES TO BE USED =',I5) 
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350 FORMAT { //5X, ' CHORDWISE BENDING STIFFNESS (LB-IN**2) ' ) 

352 FORMAT (//5X, ' TORSIONAL STIFFNESS (LB-IN**2) ' ) 

354 FORMAT (/ /5X, ' DISTANCE BETWEEN MASS AND ELASTIC AXIS(IN)=') 

356 FORMAT (//5X, 'TWIST OF THE BLADE' , /5X, ' (DEGREES) =' ) 

358 FORMAT (//5X, 'MASS MOMENT OF INERTIA ABOUT THE CHORD (LB-SEC* * 2 )' ) 

360 FORMAT (//5X, 'MASS MOMENT OF INERTIA ABOUT AN AXIS', 

$/5X, 'PERPINDICULAR TO THE CHORD THROUGH THE CENTER OF', 

$/5X, 'GRAVITY (LB-SEC**2) ' ) 

361 FORMAT (//5X, 'THE COEFFICIENTS OF THE CLEVIS COMPATIBILITY 
$FUNCTIONS ARE' ) 

362 FORMAT (//5X, 'COEFFICENTS OF THE FIRST LOAD PATH') 

363 FORMAT (//5X, 'COEFFICENTS OF THE SECOND LOAD PATH') 

364 FORMAT (//5X, 'FLAP WISE BENDING') 

365 FORMAT (//5X, 'CHORDWISE BENDING') 

366 FORMAT (//5X, 'MODE NUMBER') 

367 FORMAT (//5X, 'TWISTING' ) 

368 FORMAT (//5X, ' STRETCHING' ) 

369 FORMAT (//5X, 'EIGENVALUES ARE AS FOLLOWS') 

370 FORMAT (//5X, 'NATURAL FREQUENCIES ARE AS FOLLOWS') 

371 FORMAT (//5X, ' EIGENVECTORS ARE AS FOLLOWS') 

372 FORMAT (//5X, 'NORMALIZED EIGENVECTORS') 

C 

600 STOP 
END 

Q *★****************★*********★*★***★*★*★**★*★★*****★★**★*★*★★*★****★**★ 

SUBROUTINE EECM (A, LAMDA, VECT, HL, H, CNT, COLM, INTH) 


C 
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C 

C COMPUTES EIGENVALUES AND EIGENVECTORS OF A COMPLEX MATRIX 

C 

C 

COMPLEX A (20, 20) , LAMDA (20) , SHIFT (3) , TEMP, TEMPI, TEMP2, SINN, COSS, 
$HL(60, 60) , CNT (20) ,VECT(20,20) ,COLM(20) ,H(60, 60) 

LOGICAL TWICE, CON 

INTEGER INTH (20,2) ,R,RP1,RP2 

N=INTH (1,1) 

M=N 

NV=INTH (2,1) 

DO 1 1=1, M 
DO 1 J=1,M 

1 HL ( J, I)=A(J, I) 

NCAL=N 

SHIFT(1)=(0. 0,0.0) 

ICOUNT=0 
IF (N-2) 2,37,3 

2 LAMDA ( 1 ) =A (1,1) 

VECT (1,1) = (1. 0,1.0) 

GO TO 63 

C 

C REDUCE MATRIX TO HESSENBURG GORM 
C 

3 NM2=N-2 

DO 14 R=1 , NM2 
RP1=R+1 


RP2=R+2 
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ABIG=0 . 0 
INTH (R, 1 ) =RP1 
DO 4 I=RP1,N 

ABSSQ=REAL (HL ( I , R) ) * * 2 +AIMAG ( HL ( I , R) ) **2 
IF (ABSSQ.LE.ABIG) GO TO 4 
INTH (R, 1 ) =1 
ABIG=ABSSQ 

4 CONTINUE 
INTER=INTH(R,1) 

IF (INTER. EQ. 0.0) GO TO 14 
IF ( INTER. EQ.RP1) GO TO 7 
DO 5 I=R, N 
TEMP=HL(RP1, I) 

HL (RP1 , I ) =HL ( INTER, I ) 

5 HL ( INTER, I ) =TEMP 
DO 6 1=1, N 
TEMP=HL ( I , RPl ) 

HL ( I , RP 1 ) =HL ( I , INTER) 

6 HL ( I , INTER) =TEMP 

7 DO 8 I=RP2,N 
COLM(I)=HL(I,R) /HL (RPl , R) 

8 HL ( I , R) =COLM ( I ) 

DO 10 1=1, RPl 
TEMP= (0 . 0, 0 . 0) 

DO 9 J=RP2 , N 

9 TEMP=TEMP+HL(I, J) *COLM(J) 

10 HL ( I , RPl ) =HL ( I , RPl ) +TEMP 

DO 12 I=RP2,N 
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TEMP=(0. 0, 0. 0) 

DO 11 J=RP2,N 

11 TEMP -TEMP +HL (I, J) *COLM(J) 

12 HL (I, RP1) =HL (I, RP1) +TEMP-COLM ( I ) *HL(RP1,RP1) 
DO 13 I=RP2,N 

DO 13 J=RP2 , N 

13 HL(I, J)=HL(I, J)-COLM(I) *HL(RP1, J) 

14 CONTINUE 
C 

C CALCULATE EPSILON 
C 

EPS=0 . 0 
DO 15 1=1, N 

15 EPS=EPS+CABS (HL (1,1) ) 

DO 17 1=2, N 

SUM=0 . 0 
IM1=I-1 
DO 16 J=IM1,N 

16 SUM=SUM+CABS (HL (I, J) ) 

IF (SUM.GT.EPS) EPS=SUM 

17 CONTINUE 
FN=N 

EPS=SQRT (FN) *EPS*1 . OE-12 
IF (EPS . LT . 1 . OE-12 ) EPS-1. OE-12 
DO 18 1=1, N 
DO 18 J=1,N 

18 H ( J, I) =HL ( J, I) 

19 IF (N-2) 20,37,21 
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20 LAMDA (M) =HL ( 1 / 1 ) +SHIFT ( 1 ) 

GO TO 39 

21 MN1=M-N+1 
CON= . FALSE . 

IF (REAL (HL (N, N) ) . NE .0.0. OR. AIMAG (HL (N, N) ) . NE . 0 . 0 ) CON=.TRUE. 

IF (.NOT. CON) GO TO 22 
HL(N,N)=HL(N,N) 

IF (ABS (REAL (HL (N, N-l ) /HL (N, N) ) ) +ABS (AIMAG (HL (N, N-l ) /HL 
$ (N,N) ) ) -1.0E-12) 23,23,22 

22 IF (ABS (REAL (HL (N, N-l) ) ) +ABS (AIMAG (HL (N, N-l ) ) ) .GE.EPS) GO TO 24 

23 LAMDA (MN1 ) -HL (N, N) +SHIFT (1) 

ICOUNT=0 

N-N-l 
GO TO 19 
C 

C DETERMINE SHIFT 
C 

24 TEMP=HL (N, N) *HL (N-l , N-l ) -HL (N, N-l ) *HL (N-l , N) 

TEMP1=HL (N-l, N-l) +HL (N, N) 

SHIFT (2 ) = (TEMP1+CSQRT ( (TEMPI) **2-4 . 0* (TEMP) ) ) /2 . 0 
IF ( SHIFT ( 2 ) . NE . ( 0 . 0 , 0 . 0 ) ) GO TO 25 
SHIFT (3) =TEMP1 
GO TO 26 

25 SHIFT (3) =TEMP/SHIFT (2) 

26 IF (CABS (SHIFT (2) -HL(N,N) ) . LT . CABS (SHIFT (3 ) -HL (N, N) ) ) GO TO 27 
INDEX=3 

GO TO 28 

27 INDEX=2 
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28 IF (CABS (HL(N-l,N-2) ) . GE . EPS) GO TO 29 
LAMDA (MN1 ) =SHIFT ( 2 ) +SHIFT ( 1 ) 

LAMDA (MN1+1 ) =SHIFT (3) + SHI FT ( 1 ) 

ICOUNT=0 
N=N-2 
GO TO 19 

29 SHIFT ( 1 ) =SHIFT ( 1 ) + SHIFT ( INDEX) 

C 

C PERFORM GIVEN ROTATIONS, QR ITERATIONS 
C 

DO 30 1=1, N 

30 HL (1,1) =HL (1,1) -SHIFT ( INDEX) 

IF ( I COUNT . LE .40) GO TO 31 
NCAL=M-N 

GO TO 39 

31 NM1=N-1 
TEMP1=HL (1,1) 

TEMP 2 =HL (2,1) 

DO 36 R=1 , NM1 
RP 1 =R+ 1 

SUM=SQRT (REAL (TEMPI) **2+AIMAG (TEMPI) **2+REAL (TEMP 2 ) **2+AIMAG 
$ (TEMP2) **2 ) 

IF (SUM. NE. 0.0) GO TO 32 
TEMP1=HL (RP1 , RP1 ) 

TEMP2=HL (R+2 , RP1 ) 

GO TO 36 

32 COSS=TEMPl/SUM 
SINN=TEMP2/SUM 



INDEX=MAXO (R-l, 1) 

DO 33 I=INDEX, N 

TEMP=CONJG (COSS) *HL (R, I ) +CONJG (SINN) *HL(RP1,I) 
HL(RP1, I) — SINN*HL (R, I ) +COSS*HL (RP1 , I ) 

33 HL (R, I ) =TEMP 
TEMP 1 =HL ( RP 1 , RP 1 ) 

TEMP2=HL ( R+ 2 , RP 1 ) 

DO 34 1=1, R 

TEMP=COSS*HL ( I , R) +SINN*HL ( I , RP1 ) 

HL ( I , RP 1 ) =-CON JG (SINN) *HL (I , R) +CONJG ( COSS ) *HL ( I , RP 1 ) 

34 HL ( I , R) =TEMP 
INDEX=MINO (R+2 , N) 

DO 35 I=RP1, INDEX 

HL ( I , R) =SINN*HL ( I , RP1 ) 

35 HL ( I , RP1 ) =CONJG (COSS) *HL ( I , RP1 ) 

36 CONTINUE 
ICOUNT=ICOUNT+l 
GO TO 21 

37 TEMP1=HL (N-l , N-l ) +HL (N, N) 

TEMP2=HL(N,N) *HL (N-l , N-l ) -HL (N, N-l ) *HL(N-1,N) 

TEMP= ( TEMP 1 +CSQRT ( (TEMPI ) **2-4 . 0* (TEMP2 ) ) ) /2 . 0 
IF (REAL (TEMP ) .NE. 0.0) GO TO 38 
LAMDA(M) =SHIFT (1) 

LAMDA (M-l ) “TEMP1+SHIFT (1 ) 

GO TO 39 

38 LAMDA (M) “TEMP+SHIFT ( 1 ) 

LAMDA (M-l) “TEMP2/TEMP+SHIFT (1) 
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VALUSE IN DECENDING ORDER OF ABSOLUTE MAGITUDE 

39 CONTINUE 

IF (M. NE . 2 ) GO TO 40 

EPS=AMAX1 (CABS (LAMDA(l) ) , CABS ( LAMDA ( 2 ) ) ) *1 . OE-12 
IF (EPS. LT. 1 . OE-12) EPS-1. OE-12 
H(1,1)-HL(1,1) 

H (1, 2) =HL (1,2) 

H (2,1) =HL (2,1) 

H (2, 2) =HL (2,2) 

40 CONTINUE 

IF (NCAL.LE.l) GO TO 143 

K=NCAL+1 

DO 41 1=1 , NCAL 

K-K-l 

CNT (I) -LAMDA (K) 

41 CONTINUE 
L-NCAL-1 

DO 43 1=1, L 
K-I + l 

DO 42 J=K, NCAL 

IF (CABS (CNT (I) ) .GE.CABS (CNT (J) ) ) GO TO 42 
TEMP=CNT ( J) 

CNT ( J) =CNT (I) 

CNT (I) -TEMP 

42 CONTINUE 

43 LAMDA ( I ) =CNT ( I ) 

LAMDA (NCAL) =CNT (NCAL) 
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C 

C CALCULATES VECTORS 
C 

143 CONTINUE 

NC=MIN0 (NV,NCAL) 

IF (NC.EQ. 0) GO TO 63 

NM1-M-1 

DO 62 L=1,NC 

DO 45 1=1, M 

DO 44 J=1,M 

44 HL ( J, I ) =H ( J, I ) 

45 HL (I, I) =HL (I, I) -LAMDA(L) 

DO 49 1=1, NM1 

COLM (I) = (0. 0,0.0) 

INTH (I , 2 ) =0 
IP1-I+1 

IF (CABS (HL (1+1,1) ) . LE . CABS (HL (1,1) ) ) GO TO 47 
INTH (I, 2) =1 
DO 46 J=1,M 
TEMP-HL (IP1, J) 

HL (IP1 , J) =HL ( I , J) 

46 HL (I, J) =TEMP 

47 IF (HL (I, I) . EQ. (0 . 0, 0 . 0) ) GO TO 49 
COLM (I ) =-HL (IP1 , I ) /HL (1,1) 

DO 48 J=IP1 , M 

48 HL (IP1, J) =HL (IP1, J) +COLM(I) *HL (I, J) 

49 CONTINUE 


DO 50 1=1, M 
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50 CNT (I) = (1 . 0, 0 . 0) 

TWICE=. FALSE. 

51 IF (HL(M,M) .EQ. (0.0, 0.0) ) HL(M,M)=EPS 
CNT (M) =CNT (M) /HL (M, M) 

DO 53 1=1 , NM1 
K=M-1 

DO 52 J=K,NM1 

52 CNT (K) =CNT (K) -HL (K, J+l ) *CNT(J+1) 

IF (HL (K,K) .EQ. (0.0, 0.0)) HL(K,K)=EPS 

53 CNT (K) =CNT (K) /HL (K, K) 

ABIG=0.0 

DO 54 1=1, N 

SUM=ABS (REAL (CNT (I) ) ) +ABS (AIMAG (CNT (I) ) ) 
IF (SUM. GT . ABIG) ABIG=SUM 

54 CONTINUE 

DO 55 1=1, M 

55 CNT (I) =CNT (I) /ABIG 
IF (TWICE) GO TO 57 
DO 56 1=1, NM1 

IF ( INTH (1,2) . EQ . 0 ) GO TO 56 
TEMP=CNT (I) 

CNT ( I ) =CNT ( I + 1 ) 

CNT ( I + 1 ) =TEMP 

56 CNT (1+1) =CNT (1+1) +COLM(I) *CNT (I) 

TWICE= . TRUE . 

GO TO 51 

57 IF (M. EQ . 2 ) GO TO 60 


NM2=M-2 
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DO 59 1=1, NM2 
N1I-M-1-I 
NI1-M-I+1 
DO 58 J=NI1,M 

58 CNT ( J) =H ( J, Nil ) *CNT (N1I+1) +CNT ( J) 

INDEX=INTH (Nil, 1) 

TEMP=CNT (N1I+1) 

CNT (N1I+1) =CNT (INDEX) 

59 CNT (INDEX) =TEMP 
C 

C NORMALIZE VECTORS TO UNIT LENGTH 
C 

60 CONTINUE 
TEMP= (0 . 0, 0 . 0) 

DO 61 1=1, M 

61 TEMP=TEMP+CNT ( I ) *CNT ( I ) 

TEMP=CSQRT (TEMP) 

IF (TEMP. EQ. (0.0, 0.0) ) TEMP= ( 1 . 0 , 0 . 0 ) 

DO 62 1=1, M 

VECT ( I , L) =CNT ( I ) /TEMP 

62 CONTINUE 

63 CONTINUE 

INTH (1,1) =NCAL 

RETURN 

END 

q ******************************************************************* 

SUBROUTINE GFCTS (M) 


C 
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C 

C THIS SUBROUTINE CALCULATES THE GALERKIN'S FUNCTIONS 

C 

C 

REAL MYL (5,2,11) , MYLP (5,2,11) ,MZL(5,2,11) , MZLP (5 , 2 , 11 ) ,MXL(5,2,11) 
1 , MXLP (5,2,11) ,MY (5,51) , MYP (5, 51), MZ (5, 51) ,MZP ( 5 , 51 ) , MX (5 , 51 ) , 

2MXP (5,51) 

DIMENSION AL (5) , WL (5, 2 , 11 ) , WLP (5,2, 11)-, VL (5, 2, 11) ,VLP (5,2,11), 

1TL (5,2,11) , VXL (5,2,11) ,VXLP (5, 2, 11) , W (5, 51 ) , WP (5, 51) , V (5, 51) , 

2VP (5,51) ,T (5,51) ,U (5,51) ,VX (5,51) ,VXP (5,51) ,C (4,4) ,CT (4,4) , 
3B(5,2,4) ,CT1 (4,4) ,A(5,4,4) , FD (2,2) , ul (5,2,11) 

DATA AL/ 1 .875,4. 694,7 . 855,10. 996,14.137/ 

COMMON /XM1 /NPATH 

COMMON /XM2 /SPAN1 , SPAN 2 , SPAN3 , FD , OMEGAN 
COMMON/XGFCT1 /A, B 

COMMON/XGFCT2 /MYL , MYLP , MZ L , MZLP , MXL, MXLP , MY , MYP , MZ , MZP , MX , MXP , WL , 
1WLP , VL, VLP , TL, UL, VXL, VXLP , W, WP , V, VP , T, U, VX, VXP 
DO 220 1=1, M 
X=0 . 0 

SN=SIN (AL (I) ) 

SH=SINH (AL ( I ) ) 

CN=COS (AL (I) ) 

CH=COSH ( AL ( I ) ) 

E= (-SN+SH) / 2 . 0*SN*SH 

G= (CN+CH) / 2 . 0*SN*SH 

D= (2 . 0*FLOAT (I) -1 .0) *Pl/2 . 0 

DO 210 J-1,11 

X1=X*D 
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X=X*AL ( I ) 

SN=SIN (X) 

S H=SINH (X) 

CN=COS (X) 

CH=COSH (X) 

CN1=C0S (XI) 

SNl-SIN(Xl) 

DO 200 K=1 , NPATH 

WL(I,K, J)=E* (SN-SH) +G* (CN-CH) 

WLP ( I , K, J) =AL (I) *E* (CN-CH) +AL (I) *G* (-SN-SH) 

MYL ( I , K, J) =AL ( I ) **2*E* (-SN-SH) +AL (I) **2*G* (-CN-CH) 
MYLP(I,K, J)-AL(I) **3*E* ( -CN-CH) +AL ( I ) * * 3*G* (SN-SH) 
VL ( I , K, J) =WL ( I , K, J) 

VLP (I , K, J) =WLP ( I , K, J) 

MZL(I,K, J) =MYL (I,K, J) 

MZLP ( I , K, J) =MYLP ( I , K, J) 

TL ( I , K, J) =SN1 
MXL (I , K, J) =D*CN1 
MXLP (I,K, J) =-D*D*SNl 
UL (I, K, J) =TL (I, K, J) 

VXL ( I , K, J) =MXL ( I , K, J) 

VXLP ( I , K, J) =MXLP ( I , K, J) 

200 CONTINUE 
X=X+HS1 
210 CONTINUE 
X=SPAN1 
DO 220 J-1,51 


X1=X*D 
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X=AL ( I ) 

SN=SIN (X) 

SH=SINH (X) 

CN=COS (X) 

CH=COSH (X) 

CN1=C0S (XI) 

SN1=SIN (XI) 

W(I, J)-E* (SN-SH) +G* (CN-CH) 

WP (I, J)-AL(I) *E* (CN-CH) +AL (I) *G* (-SN-SH) 

MY (I, J) =AL (I) **2*E* (-SN-SH) +AL ( I ) **2*G* (-CN-CH) 
MYP (I, J)-AL(I) **3*E* (-CN-CH) +AL(I) **3*G* (SN-SH) 
V (I , J) =W (I , J) 

VP (I, J) =WP (I, J) 

MZ (I, J) =MYP (I, J) 

T(I, J)-SN1 
MX(I, J)=D*CN1 
MXP (I, J) =-D*D*SNl 
U (I, J) — T (I, J) 

VX (I, J) =MX (I, J) 

VXP (I, J) =MXP (I, J) 

X=X+HS2 
220 CONTINUE 
C( 1,1) -1.0 
C (1, 2) =SPAN3 
C (1, 3) =SPAN3*SPAN3 
C (1, 4) =C (1,3) *SPAN3 
C (2 , 1) =0 . 0 
C (2, 2) *=1 . 0 
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C (2, 3) =2 . 0*SPAN3 

C (2, 4) =3. 0*SPAN3*SPAN3 

C (3,1) -1.0 

C(3,2)=SPAN1 

C (3,3) =SPAN1 * SPAN1 

C(3,4)=C(3,3)*SPAN1 

C (4, 1) =0 . 0 

C (4 , 2) =1 . 0 

C (4, 3) =2 . 0*SPAN1 

C (4, 4) =3. 0*SPAN1*SPAN1 

CALL SOLUTN(C,4,-l,4) 

N=NPATH+1 
n2=2*npath 
DO 240 1=1, M 
DO 230 K=1 , N2 , 2 
Kl-K 

IF(K.GT.l) Kl-K-1 
CT ( 1 , K ) =WL ( I , K1 , 10) 

CT (2 , K) =WLP (I,K1, 10) 

CT (3,K) =W (1,1) +FD (Kl, 1) *T (I, 1) 

CT (4 , K) =WP (1,1) 

CT (1 , K+l ) =VL (I , Kl , 10) 

CT (2, K+l) =VLP (I , Kl , 10) 

CT ( 3 , K+ 1 ) =V ( 1 , 1 ) -FD (Kl, 2) *T (I, 1) 

CT (4 , K+l ) =VP (1,1) 

B ( I , 1 , K) = (SPAN1*TL ( I , Kl , 10 ) -SPAN3*T (1,1) ) / ( SPAN1-SPAN3 ) 
B(I,2,K)=(T(I,1) -TL ( I , Kl , 10 ) ) / ( SPAN1 -SPAN3 ) 

B ( I , 1 , K+l ) = ( SPAN1 *UL ( I , Kl ,10) -SPAN3*U (1,1) ) / (SPAN1-SPAN3) 
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B ( I , 2 , K+l ) = (U ( I , 1 ) -UL ( I , K1 , 10 ) ) / (SPAN1-SPAN3 ) 

230 CONTINUE 

CALL MATMUL(C,CT,CT1,4,4,N2) 

DO 240 11=1,4 
DO 240 I J=1 , N2 
A(I,II,IJ) =CT1 (II, IJ) 

240 CONTINUE 

DO 260 1=1, MODES 
DO 260 J=1 , NPATH 
I J=2* J-l 
TL ( I , J, 11 ) =0 . 0 
UL ( I , J, 11 ) =0 . 0 
WL (I, J, 11) =0 . 0 
VL(I, J, 11) =0.0 
WLP (I, J, 11) =0.0 
VLP(I, J, 11) =0.0 

WL ( I , J , 1 1 ) =A ( I , L , I J ) *SPAN1** ( L- 1 ) +WL ( I , J, 1 1 ) 

VL ( I , J, 11 ) =A ( I , L, I J+l ) *SPAN1** (L-1)+VL(I, J, 11) 

WLP (I , J, 11 ) =FLOAT (L-l ) *A ( I , L, I J) *SPAN1** (L-2J+WLP (I , J, 11 ) 
VLP ( I , J, 11 ) =FLOAT (L-l ) *A ( I , L, I J+l ) *SPAN1** (L-2)+VLP (I , J, 11 ) 
250 CONTINUE 

TL(I, J, 11) =B (I, 1, IJ) +B (I, 2, IJ) *SPAN1 
UL (I, J, 11) =B (I, 1, I J+l) +B (I, 2, I J+l) *SPAN1 
260 CONTINUE 
RETURN 
END 


SUBROUTINE GPROC (MODES) 
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THIS SUBROUTINE FORMULATES THE EIGENVALUE PROBLEM 


REAL INTG 

DIMENSION All (5, 5) ,A12 (5,5),A21(5,5) ,A22 (5, 5), A3 1 (5, 5) ,A32 (5,5) , 
1A33 (5,5) , A41 (5,5) ,A13 (5,5) ,A14 (5,5) ,A23 (5,5) ,A24 (5,5) ,A34 (5,5) , 

2 A3 5 (5,5) , A36 (5,5) ,A42 (5, 5) , Bll (5, 5) , B12 (5, 5) , B21 (5, 5) ,b33 (5,5) , 
3B22 (5,5) , B31 (5,5) ,B32 (5,5),B41(5,5) ,A51 (5,5) ,A52 (5,5) ,A61 (5,5) , 
4A62 (5,5) ,A71 (5,5) ,A72 (5,5) ,A73 (5,5) ,A83 (5 , 5 ) , Dll ( 2 , 11 ) , D12 (2 , 11 ) , 
5D13 (2,11) , D14 (2,11) ,D15 (2, 11) , D16 (2, 11) , D17 (2, 11) ,D18 (2,11), 

6D19 (2,11) , DUO (2, 11) ,D21 (51) ,D22 (51) ,D23 (51) ,dlll (2, 11) , 

7D24 (51) , D25 (51) ,D26 (51) ,D27 (51) ,D28 (51) ,D2 9 (51) ,D210 (51) , D211 (51) , 
8D212 (51) ,F11 (2, 11) ,F12 (2,11) ,F13 (2 , 11 ) , F14 (2 , 11 ) ,F15 (2,11) ,F21 
9(51) ,F22 (51) ,F23(51) ,F24 (51) ,F25 (51) , FD (2 , 2 ) , F16 (2 , 11 ) , F17 (2,11) , 
1F26 (51) ,F27 (51) ,D1 (5,5) ,D2 (5,5) ,C(10,10) , Cl (10, 10) ,C2 (10,10) , 
2WL(5,2, 11) , WLP (5,2, 11) ,VL (5, 2, 11) ,VLP (5, 2 , 11 ) , TL (5 , 2, 11 ) , UL (5, 2 , 
311) ,VXL(5,2, 11) , VXLP (5, 2, 11) , W (5, 51) , WP (5, 51) , V (5, 51) , VP (5, 51) , 
4t(5,51) , U ( 5 , 51 ) , VX (5,51) ,VXP (5,51),A(5,4,4),B(5,2,4),E1(20,20), 

5e2 (20,20) ,E3 (20,20) ,E4 (20,20) ,DS (20, 20) ,dll2(2,ll) , dll3 (2, 11) , 
6d213 (51) , a81 (5, 5) ,a82 (5, 5) 

REAL MYL (5,2,11) , MYLP (5,2,11) ,MZL(5,2,11) ,MZLP (5,2,11) , 

Imxl (5,2,11) , MXLP (5,2,11) ,MY(5,51) ,MYP (5, 51) ,MZ (5, 51) ,MZP(5,51) , 

2mx (5, 51) , MXP (5,51) 

COMMON /XM1 /NPATH 

COMMON/XM2 / SPAN1 , SPAN 2 , SPAN3 , FD , OMEGAN 

COMMON / XM3 /Dll, D12,D13,D14,D15,D16,D17,D18,D19,D110, Dill, D112,D113 
1,D21,D22,D23,D24,D25,D26,D27,D28,D29,D210,D211,D212,D213 
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C0MM0N/XGFCT2 /MYL , MYLP , MZL, MZLP , MXL, MXLP , MY , MYP , MZ , M2P , MX , MXP , WL , 
1WLP,VL,VLP,TL,UL,VXL,VXLP,W,WP,V,VP,T,U,VX,VXP 
C0MM0N/XGPR1/DS 
DO 120 K=l, MODES 
DO 120 J=l, MODES 
DO 100 1=1 , NPATH 
DO 100 L-1,11 

Fll ( I , L) =D113 (I,L) *WLP ( J, I , L) *WLP(K,I,L) 

F12 (I,L)=D12 (I,L) *TL(J, I,L) *WLP (K,I,L) 

F13 (I, L)— MYLP ( J, I , L) *WLP (K,I,L) 

F14 (I,L) -Dll (I, L) *WL ( J, I,L) *WL(K, I,L) 

F15 ( I , L ) =D 1 3 (I,L) *TL(J,I, L) *WL(K,I,L) 

100 CONTINUE 

DO 110 L-1,51 

F21 (L)=D213 (L) *WP (J,L) *WP (K, L) 

F22 (L) =D22 (L) *T ( J, L) *WP (K, L) 

F23 (L) =-MYP ( J, L) *WP (K, L) 

F24 (L)=D21 (L) *W(J,L) *W(K,L) 

F25 (L)=D23 (L) *T(J,L) *W(K,L) 

110 CONTINUE 

All (K, J) =INTG (Fll , FI 2 ) 

A12 (K, J)=INTG(F12,F22) 

A14 (K, J)=INTG(F13,F23) 

Bll (K, J)=INTG(F14,F24) 

B12 (K, J)=INTG(F15,F25) 

A14 (K, J) =A14 (K, J) +W (K, 1 ) * (MYL ( J, 1, 11) +MYL ( J, 2, 11) -MY ( J, 1) ) 

A13 (K, J) =-FD (1,2) *VXL(J, 1, 11) -FD (2, 2) *VXL (J, 2, 11) 


120 CONTINUE 
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DO 150 K=l, MODES 
DO 150 J-l, MODES 
DO 130 1=1 , NPATH 
DO 130 L=1 ,11 

Fll (I , L) =D113 ( I , L) *VLP (J, I,L) *VLP (K, I , L) 

F12 (I , L) =D14 ( I , L) *TL(J, I,L) *VLP (K, I , L) 

F13 (I,L) =-MZLP ( J, I, L) *VLP (K,I,L) 

F14 (I,L) -Dll (I,L) *VL(J, I,L) *VL (K, I, L) 

F15 (I,L)— D15 (I,L) *TL ( J, I, L) *VL (K, I , L) 

130 CONTINUE 

DO 140 L-1,51 

F21 (L)=D213 (L) *VP (J,L) *VP (K, L) 

F22 (L)=D24 (L) *T ( J,L) *VP (K, L) 

F23 (L) =-MZP ( J, L) *VP (K, L) 

F24 (L)=D21 (L) *V(J,L) *V(K,L) 

F25 (L) =-D25 (L) *T ( J, L) *V(K,L) 

140 CONTINUE 

A21 (K, J) =INTG (Fll , F21 ) 

A22 (K, J) =INTG (F12 , F22 ) 

A24 (K, J) =INTG (F13,F23) 

B21 (K, J) =INTG (F14 , F24 ) 

B22 (K, J) =INTG (F15, F25) 

A24 (K, J) =A24 (K, J) +V(K, 1) * (MZL ( J, 1 , 11 ) +MZL ( J, 2 , 11 ) -MZ ( J, 1) ) 
A23 (K, J) =-FD (1,1) *VXL(J, 1, 11) -FD (2 , 1 ) *VXL ( J, 2 , 11 ) 

150 CONTINUE 

DO 180 K=l, MODES 
DO 180 J-l, MODES 
DO 160 1=1, NPATH 


C 
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DO 160 L-1,11 

Fll (I, L) — D12 ( I , L) *VLP ( J, I,L) *TL (K, I,L) +OMEGAN*D15 (I,L) * 
1VL ( J, I , L) *TL (K, I , L) 

F12 ( I , L) =D14 (I, L) *WLP (J,I,L) *TL (K, I , L) 

F13 (I,L)-D16(I,L) *TL(J, I , L) *TL (K, I , L) 

F14 (I,L)=-MXLP ( J, I , L) *TL (K, I , L) 

F15 (I,L) -D17 (I,L) *TL ( J, I , L) *TL(K, I,L) 

FI 6 (I,L)-D13 ( I , L) *WL(J, I,L) *TL(K,I,L) 

F17 ( I , L) =-D15 (I, L) *VL ( J, I , L) *TL(K, I,L) 

160 CONTINUE 

DO 170 L-1,51 

F21 (L) =-D22 (L) *VP ( J,L) *T (K, L) +OMEGAN*D25 (L) *V( J,L) *T (K, L) 
F22 (L)=D24 (L) *WP (J,L) *T(K,L) 

F23 (L)=D2 6(L) *T(J,L) *T(K,L) 

F24 (L)=-MXP ( J, L) *T(K,L) 

F25 (L) =D27 (L) *T ( J,L) *T (K,L) 

F2 6 (L)=D23 (L) *W( J,L) *T (K,L) 

F27 (L)«- D25 (L) *V(J,L) *T(K,L) 

170 CONTINUE 

A31 (K, J) =INTG (Fll , F12 ) 

A3 2 (K, J) =INTG (F12 , F22 ) 

A3 3 (K, J) =INTG (F13,F23) 

A3 6 (K, J)=INTG(F14,F24) 

B31 (K, J) =INTG (F17 , F27 ) 

B32 (K, J) =INTG (F16,F26) 

B33 (K, J) =INTG (F15 , F25 ) 

A3 6 (K, J) =A36 (K, J) +T (K, 1 ) * (MXL ( J, 1, 11) +MXL ( J, 2, 11) -MX ( J, 1) ) 
A3 5 (K, J) =A35 (K, J) -FD (1, 1) *MYLP ( J, 1, 11) *T (K, 1) -FD (2,1) 
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1 *MYLP ( J, 2,11) *T(K,1) 

A3 2 (K, J) =A32 (K, J) +T(K, 1) * (D113 (1, 11) *WLP (J, 1,11) *FD (1,1) 
1+D113 (2,11) *WLP ( J, 2, 11) *FD(2,1) ) 

A3 3 (K, J) =A33 (K, J) +T (K, 1) * (FD (1,1) *D12 (1, 11) *T(J, 1) ) 

1+FD (2, 1) *D12 (2,11) *T(J,2)+T(K,1) * (FD(2,1) *D14 (2,11) * 

2TL ( J, 1,11) +FD (2,2) *D14 (2, 11) *TL ( J, 2, 11) ) 

A3 4 (K, J ) =A3 4 (K, J) +T (K, 1 ) * (FD ( 1 , 2 ) *MZLP ( J, 1, 11) +FD (2,2) 
1*MZLP ( J, 2,11) ) 

A31 (K, J) =A31 (K, J) *T (K, 1) * (-FD (1, 2) *D113 (1, 11) *VLP (J, 1, 11) 
1-FD (2, 2) *D113 (2,11) *VLP ( J, 2, 11) ) 

180 CONTINUE 

DO 210 K— 1 , MODES 
DO 210 J-l, MODES 
DO 190 1=1 , NPATH 
DO 190 L-1,11 

Fll ( I , L) =-OMEGAN*Dl 1 (I , J) *UL ( J, I, L) *UL(K,I,L) 

F12 ( I , L) =-VXLP ( J, I , L) *UL(K, I,L) 

F13 (I , L) =D11 ( I , L) *UL ( J, I , L) *UL(K, I,L) 

190 CONTINUE 

DO 200 L-1,51 

F22 (L) =-OMEGAN*D21 (L) *U ( J, L) *U (K, L) 

F22 (L)— VXP ( J, L) *U(K,L) 

F23 (L)=D21 (L) *U ( J, L) *U (K, J) 

200 CONTINUE 

A41 (K, J) =INTG (Fll , F21 ) 

A42 (K, J) =INTG (F12 , F22 ) 

B41 (K, J) =INTG (F13 , F23 ) 

A42 (K, J) =A42 (K, J) +U (K, 1) * (VXL ( J, 1, 11) +VXL ( J, 2, 11) -VX ( J, 1) ) 



210 CONTINUE 


DO 240 K=l, MODES 
DO 240 J=1 , MODES 
DO 220 1=1 , NPATH 
DO 220 L-1,11 

Fll (I, L) =UL ( J, I , L) *VXL (K, I , L) 

F12 (I,L)=-D18 (I, J) *VXL ( J, I , L) *VXL(K, I,L) 
F13 (I,L)-TL(J,I,L) *MXL(K / I / L) 

F14 (I, L) =-D19 (I, L) *MXL(J,I,L) *MXL(K, I,L) 
220 CONTINUE 

DO 230 L-1,51 

F21 (L)=U(J,L) *VX(K,L) 

F22 (L)=-D28 (L) *VX(J,L) *VX(K,L) 

F23 (L)=T ( J, L) *MX (K, L) 

F24 (L) =-D2 9 (L) *MX(J,L) *MX(K,L) 

230 CONTINUE 

A51 (K, J) =INTG (Fll , F21 ) 

A52 (K, J)=INTG(F12,F22) 

A61 (K, J) =INTG (F13,F23) 

A62 (K, J)=INTG(F14,F24) 

240 continue 

DO 270 K=l, MODES 
DO 270 J-l, MODES 
DO 250 1=1, NPATH 
DO 250 L-1,11 

Fll (I, L) =MZL ( J, I, L) *MZL(K, I,L) 

F12 (I, L) -Dill (I, L) *MYL ( J, I , L) *MZL(K, I,L) 
F13 (I,L)=- D110 ( I , L) *MYL ( J, I , L) *MZL(K, I,L) 
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F14 (I,L)=MYL(J, I, L) *MYL (K, I , L) 

F15 (I,L)=-D112 ( I , L) *MYL ( J, I , L) *MYL (K, I , L) 
FI 6 (I,L) =D112 ( I , L) *MZL ( J, I, L) *MYL (K, I , L) 
250 CONTINUE 

DO 260 L=1 , 51 

F21 (L)=MZ ( J, L) *MZ (K, L) 

F22 (L) =D211 (L) *MY(J,L) *MZ (K, L) 

F23 (L) =-D210 (L) *MY(J,L) *MZ (K,L) 

F24 (L) =MY ( J, L) *MZ (K, L) 

F25 (L) =-D112 (2 , L) *MY(J,L) *MY (K, L) 

F26 (L) =D112 (2 , L) *MY (J,L) *MZ (K, L) 

260 CONTINUE 

A71 ( J, K) =INTG (Fll , F21 ) 

A72 ( J, K) =INTG (F13, F23) 

A73 ( J, K) =INTG (F12 , F22 ) 

A81 (J,K)=INTG(F14,F24) 

A82 ( J, K) =INTG (F16, F26) 

A83 ( J, K) =INTG (F15, F25) 

270 CONTINUE 

CALL SOLUTN (A52, MODES, -1, MODES) 

CALL MATMUL (A52,A51,D1, 5, 5, 5) 

CALL SOLUTN (A62, MODES, -1, MODES) 

CALL MATMUL (A62,A61, D2, 5, 5, 5) 

DO 280 1*1 , MODES 
DO 280 J-l, MODES 
C(I, J)=A72 (I, J) 

C (I, J+MODES) =A73 (I, J) 

C (I+MODES, J) =A82 (I, J) 



C (I+MODES, J+MODES) =A8 3 (I, J) 

Cl (I, J)=A71 (I, J) 

Cl (I, J+MODES) =0.0 
Cl (I+MODES, J) “0.0 
Cl (I+MODES, J+MODES) “0.0 
280 CONTINUE 
N2=2*MODES 

CALL SOLUTN (C,N2,-1,N2) 

CALL MATMUL (C, Cl , C2 , N2 , N2 , N2 ) 
N3=3*MODES 
N4=4*MODES 
DO 285 I“1,N4 
DO 285 J=1,N4 
El (I, J) =0 . 0 
E2 (I , J) =0 . 0 
E3 (I , J) =0 . 0 
E4 (I, J)=0.0 
285 CONTINUE 

DO 290 1=1, MODES 
DO 290 J=l, MODES 
El (I, J+N2) “All (I, J) 

El (I, J+N3) =A12 (I, J) 

El (I+MODES, J+MODES) =A21 (I, J) 
El (I+MODES, J+N3)=A2 2 (I, J) 

El (I+N2, J+MODES) =A31 (I, J) 

El (I+N2, J+N2) =A32 (I, J) 

El (I+N2, J+N3) =A3 3 (I, J) 

El (I+N3, J)=A41 (I, J) 
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E2 (I, J)«A13 (1/ J) 

E2 (I, J+N2)=A14 (I, J) 

E2 (I+MODES, J)=A23 (I, J) 

E2 (I+MODES, J+MODES)=A2 4 (I, J) 

E2 (I+N2, J+MODES) =A34 (I, J) 

E2 (I+N2, J+N2) =A35 (I, J) 

E2 (I+N2 , J+N3) =A36 (I, J) 

E2 (I+N3, J)=A42 (I, J) 

E3 (I, J+N2)=B11 (I, J) 

E3 (I, J+N3) =B12 (I, J) 

E3 (I+MODES, J+MODES) =B21 (I, J) 

E3 (I+MODES, J+N3)=B22 (I, J) 

E3 ( I+N2 , J+N2 ) =B32 ( I , J) 

E3 ( I +N2, J+MODES )=B31 (I, J) 

E3 (I+N2, J+N3) =B33 (I, J) 

E3 (I+N3, J)=B41 (I, J) 

E4 (I, J)— D1 (I, J) 

E4 (I+N3, J+N3)=-D2 (I, J) 

E4 (I+MODES, J+MODES) =C2 (I, J) 

E4 (I+MODES, J+N2) =C2 (I, J+MODES) 

E4 (I+N2, J+MODES) =C2 (I+MODES, J) 

E4 ( I+N2 , J+N2 ) =C2 ( I+MODES , J+MODES) 
290 CONTINUE 

CALL MATMUL (E2,E4,DS,N4,N4,N4) 
CALL SOLUTN (E3,N4,-1,N4) 

DO 300 1=1, N4 
DO 300 J=1,N4 
E2 (I, J)=E1 (I, J) +DS (I, J) 
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300 CONTINUE 

CALL MATMUL(E3,E2,DS,N4,N4,N4) 

RETURN 

END 

FUNCTION INTG (FI , F2 ) 

Q ★★★**★★★***★***★**★★★*★*★★★*★**★*******★*★★*★★******★****★*******★ 

c 

C THIS FUNCTION INTEGRATES THE FUNCTION IN THE DOMAIN 

C 

C 

REAL INTG 

DIMENSION FI (2,11) ,F2 (51) 

COMMON / XM1 / NPATH 
COMMON / XM5 / HS1,HS2 
Hl=HSl/24 . 0 
H2=HS2/24 . 0 
SUM=0 . 0 

DO 110 J-l, NPATH 

SUM=SUM+ (9. 0*F1 (J,1)+19.0*F1 {J,2)-5.0*F1 (J, 3)+Fl (J,4) ) *H1 
DO 100 K-1,8 

SUM=SUM+H1 * ( -FI ( J,K) +13. 0*F1 (J,K+1)+13. 0*F1 ( J,K+2) -FI ( J,K_3) ) 

100 CONTINUE 

SUM=SUM+H1 * (FI (J, 8) -5. 0*F1 (J, 9) +19. 0*F1 ( J, 10) +9. 0*F1 ( J, 11) ) 

110 CONTINUE 

SUM=SUM+(9.0*F2 (1) +19 . 0*F2 (2) -5 . 0*F2 (3)+F2 (4) ) *H2 
DO 120 K-1,48 

SUM=SUM+H2* (-F2 (K) +13 . 0*F2 (K+l) +13 . 0*F2 (K+2) -F2 (K+3) ) 
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120 CONTINUE 

INTG=SUM+H2* (F2 (48) -5 . 0*F2 (4 9) +19. 0*F2 (50) +9. 0*F2 (51) ) 

RETURN 

END 

SUBROUTINE INTPOL (N, A, H) 


C THIS SUBROUTINE INTERPOLATES FOR THE REQUIRED VALUES 


DIMENSION A (51) ,STA(51) , TABLE (51,1) ,B(51) 
COMMON/XM4 / STA, NS 
NN-N-1 
A (N) =A (NS ) 

NM1-NS-1 
DO 20 1*1, NM1 

20 TABLE (I,1)=(A(I+1) -A ( I ) ) / (STA ( 1+1 ) -STA ( I ) ) 
XARG=H 
DO 35 1*2, NN 
DO 25 J=1 , NS 

IF ( J . EQ . NS . OR . XARG . LE . STA ( J) ) GO TO 30 
25 CONTINUE 
30 MAX* J 

IF (MAX. LE . 2 ) MAX* 2 
ISUB=MAX-1 
YEST=TABLE (I SUB, 1) 

B ( I ) =YEST* (XARG-STA (ISUB) )+A(ISUB) 

35 XARG=XARG+H 
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DO 40 J=2,NN 
40 A(J)=B(J) 

RETURN 

END 

C ***************************************************************** 

SUBROUTINE MATMUL (A, B, C, L,M, N) 

Q ***************************************************************** 

c 

C MATRIX MULTIPLICATION 

C 

C 

DIMENSION A(L,M) ,B(M,N) ,C(L,N) 

DO 100 1=1, L 
DO 100 J=1 , N 
C(I,J)=0.0 
DO 100 K=1 , M 

C(I, J)=C(I, J) +A(I,M) *B(M, J) 

100 CONTINUE 
RETURN 
END 

Q *************************************************************** 

SUBROUTINE SOLUTN (A, N, INDIC, NRC) 

Q ★★*★★★★*★*****★***★★*★★★★★★★★★★★★★★★★★★*★★**********★***★★**★★* 

c 

C SOLUTION OF SIMULTANEOUS EQUATIONS & INVERSE 

C 

C 

DIMENSION A (NRC, NRC) ,X1 (100) , IROW(IOO) , JCOL(IOO) , JORD (100) , 



1Y (100) 

IND=0 

MAX=N 

IF (INDIC.GE.O) MAX=N+1 
DETER=1 . 0 
DO 80 K=1 , N 
KM1=K-1 
PIVOT=0 . 0 
DO 60 1=1, N 
DO 60 J=1,N 
IF (K.EQ.l) GO TO 55 
DO 50 ISCAN=1, KM1 
DO 50 JSCAN=1 , KM1 
IF (I . EQ . IROW ( I SCAN) ) GO TO 60 
IF ( J.EQ. JCOL (JSCAN) ) GO TO 60 
50 CONTINUE 

55 IF (ABS (A(I,J) ) .LE. ABS (PIVOT) ) GO TO 60 
PIVOT=A ( I , J) 

IROW (K) =1 
JCOL (K) =J 

60 CONTINUE 

IF (ABS (PIVOT) .GT.0.1E-15) GO TO 65 

DETER=0 . 0 

IND-1 

WRITE (50, 61) 

61 FORMAT (' MATRIX IS ALGORITHMICALLY SINGULAR') 
RETURN 


65 IROWK=IROW(K) 



JCOLK=JCOL(K) 

DETER=DETER*P IVOT 
DO 70 J=1 , MAX 

70 A (IROWK, J) =A (IROWK, J) /PIVOT 
A ( IROWK, JCOLK) =1 . 0 /PIVOT 
DO 80 1=1, N 
AIJCK=A(I, JCOLK) 

IF (I.EQ.IROWK) GO TO 80 
A (I, JCOLK) =-AI JCK/PIVOT 
• DO 75 J=1 , MAX 

75 IF (J.NE. JCOLK) A ( I , J) =A ( I , J) -AI JCK*A (IROWK, J) 
80 CONTINUE 
DO 85 1=1, N 
IROWI=IROW ( I ) 

JCOLI= JCOL ( I ) 

JORD (IROWI) =JCOLI 

IF (INDIC.GE.0) XI (JCOLI)=A (IROWI, MAX) 

85 CONTINUE 
INTCH=0 
NM1=N-1 
DO 90 1=1, NM1 
IP1-I+1 
DO 90 J=IP1 , N 

IF (JORD (J) .GE. JORD (I) ) GO TO 90 
JTEMP= JORD ( j ) 

JORD ( J) =JORD (I) 

JORD ( I ) = JTEMP 


INTCH=INTCH+1 
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90 CONTINUE 

IF (INTCH/2*2 .NE. INTCH) DETER=-DETER 

IF (INDIC.LE.O) GO TO 94 

RETURN 

94 DO 100 J=1 , N 
DO 95 1=1, N 
IROWI=IROW ( I ) 

JCOLI= JCOL ( I ) 

Y ( JCOLI ) -A ( I ROW I , J) 

95 CONTINUE 

do 100 1=1, N 
A(I, J)-Y(I) 

100 CONTINUE 

DO 110 1=1, N 
DO 105 J=1,N 
IROWJ=IROW(J) 

JCOLJ=JCOL(J) 

Y (IROWJ)-A(I, JCOLJ) 

105 CONTINUE 

DO 110 J=1,N 
A(I, J)-Y(J) 

110 CONTINUE 
RETURN 


END 



